31:T58a,{ const ctx = DOM.context2d(200, 200); let angle = 0; let left, prevLeft; const side = html`left` const view = html`${tex`r`} is to the ${side} of ${tex`\overrightarrow{pq}`}
${ctx.canvas}` function dot(x, y, color, label) { ctx.lineWidth = 1; ctx.fillStyle = color; ctx.beginPath(); ctx.arc(x, y, 6, 0, 2 * Math.PI); ctx.stroke(); ctx.fill(); ctx.fillStyle = 'black'; ctx.font = '16px sans-serif'; ctx.fillText(label, x + 12, y - 5); } function draw() { angle += 0.03; ctx.clearRect(0, 0, 200, 200); ctx.beginPath(); ctx.moveTo(0, 0); ctx.lineTo(200, 200); ctx.stroke(); ctx.fillStyle = 'black'; ctx.moveTo(190, 196); ctx.lineTo(200, 200); ctx.lineTo(196, 190); ctx.fill(); const rx = 100 + 30 * Math.sin(angle); const ry = 100 + 30 * Math.cos(angle); ctx.beginPath(); ctx.moveTo(40, 40); ctx.lineTo(rx, ry); ctx.setLineDash([3, 3]); ctx.strokeStyle = '#aaa'; ctx.stroke(); ctx.setLineDash([]); ctx.strokeStyle = 'black'; left = orient(rx, ry, 160, 160, 40, 40) < 0; dot(rx, ry, left ? 'red' : 'blue', 'r'); dot(40, 40, 'lightgreen', 'p'); dot(160, 160, 'lightgreen', 'q'); if (left !== prevLeft) side.innerHTML = left ? 'left' : 'right'; prevLeft = left; } return animate(view, invalidation, draw); }32:T460,md`When we use [floating point arithmetic](https://en.wikipedia.org/wiki/Floating-point_arithmetic), real numbers can't be represented exactly and are rounded to a certain precision, leading to accumulation of small errors. This causes all kinds of unexpected problems when running computational geometry algorithms on floating point input. When using the simple test in the ridiculous case above, many points that are either left or right are incorrectly classified as on the line, and vise versa. Currently the most widely used solution to this problem is [robust predicates by Jonathan Shewchuk](https://www.cs.cmu.edu/~quake/robust.html), implemented with a few thousand lines of very complicated C code way above my comprehension. I got so fascinated with this topic that I wrote a modern, fast JavaScript port, [robust-predicates](https://github.com/mourner/robust-predicates). Check it out! P.S. Jacob Rus collected a [huge gallery of beautiful plaid designs](https://observablehq.com/@jrus/non-robust-arithmetic-gallery) derived from this notebook. Extra points of awesome for printing this on a real fabric!`33:[] 2f:{"title":"Robust Arithmetic in Computational Geometry","cells":[{"id":0,"value":"md`# Robust Arithmetic in Computational Geometry\n\n**An interactive exploration of numeric robustness issues based on [Robust Arithmetic Notes](https://github.com/mikolalysenko/robust-arithmetic-notes) by [Mikola Lysenko](https://twitter.com/MikolaLysenko) + demo of [robust-predicates](https://github.com/mourner/robust-predicates).**`","pinned":false,"hidden":false,"mode":"ojs"},{"id":696,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":755,"value":"md`For a point ${tex`r`} varying across a tiny range, let's plot where it lies relative to a diagonal line — to the left, right, or on the line. We'll define the line by points ${tex`p`} and ${tex`q`}, and get orientation from the sign of this expression (a seemingly trivial computation):\n\n${tex.block`(p_x - r_x)(q_y - r_y) - (p_y - r_y)(q_x - r_x)`}\n\nPlay with the parameters to see how the plot changes:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":644,"value":"viewof z = fancyRange({\n width: 400, min: 40, max: 44, value: 42,\n label: v => html`${tex`w = 2^{${-v}}`} (zoom)`\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":146,"value":"viewof a = fancyRange({\n width: 400, min: 0, max: 3, value: 0.95,\n label: v => tex`r = (x, y), x, y \\in ${v} + [0..w)`\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":154,"value":"viewof b = fancyRange({\n width: 400, min: 0, max: 30, value: 18,\n label: v => tex`q = (k, k), k = ${v}`\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":160,"value":"viewof c = fancyRange({\n width: 400, min: 0, max: 30, value: 16.8,\n label: v => tex`p = (k, k), k = ${v}`\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":37,"value":"{\n const w = Math.min(512, width);\n const ctx = DOM.context2d(w, w, 1);\n ctx.canvas.style.imageRendering = 'pixelated';\n ctx.fillRect(0, 0, w, w);\n const imgData = ctx.getImageData(0, 0, w, w);\n const size = 2 ** -z;\n \n for (let y = 0; y < w; y++) {\n for (let x = 0; x < w; x++) {\n const rx = a + x * size / w;\n const ry = a + y * size / w;\n const o = Math.sign(chosenOrient(rx, ry, b, b, c, c));\n imgData.data[4 * (y * w + x) + o + 1] = 255;\n }\n }\n ctx.putImageData(imgData, 0, 0);\n \n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":536,"value":"md`Legend: \nto the left, \non the line, \nto the right. Beautiful! Wait, _what_?!\n\n### Choose orientation tests method:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":460,"value":"viewof chosenOrient = {\n const options = {\n 'Simple': orient,\n 'Symmetrical': orientSymmetrical,\n 'Symmetrical 2': orientSymmetrical2,\n 'More stable': orientBetter,\n 'Fully robust': orientRobust\n };\n const div = html`

`;\n for (const id in options) {\n const input = html``;\n const label = html`${input} ${id}
`;\n input.oninput = () => {\n div.value = options[id];\n div.dispatchEvent(new CustomEvent('input'));\n };\n if (id === 'Simple') input.checked = true;\n div.appendChild(label);\n }\n div.value = orient;\n return div;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1007,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":961,"value":"md`## Orientation tests explained\n\n#### Simple\nA naive implementation of the formula above.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":7,"value":"orient = (rx, ry, qx, qy, px, py) => (px - rx) * (qy - ry) - (py - ry) * (qx - rx)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1018,"value":"md`\n#### Symmetrical\n\nInstead of\n${tex.block`(p - r)\\wedge(q - r)`}\nwe can equivalently use the symmetrical expression\n${tex.block`p\\wedge q + q\\wedge r + r\\wedge p`}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1017,"value":"orientSymmetrical = (rx, ry, qx, qy, px, py) =>\n (px * qy - py * qx) + (qx * ry - qy * rx) + (rx * py - ry * px)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1016,"value":"md`\nThe below variant has much worse performance but is included because it makes quite an interesting picture. Did we mention that floating point addition is not associative?\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1015,"value":"orientSymmetrical2 = (rx, ry, qx, qy, px, py) =>\n (px * qy + qx * ry + rx * py) - (py * qx + qy * rx + ry * px)","pinned":true,"hidden":false,"mode":"ojs"},{"id":977,"value":"md`#### More stable\nA little more sophisticated test that borrows error bounds check from J. Shewchuk. If we're confident in the orientation, return it, otherwise test again with rotated point order, doing 3 combinations overall. This makes the test more precise and also symmetric — the sign will be the same regardless of the order of points.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":854,"value":"orientBetter = {\n function orient(px, py, rx, ry, qx, qy) {\n const l = (ry - py) * (qx - px);\n const r = (rx - px) * (qy - py);\n return Math.abs(l - r) >= 3.3306690738754716e-16 * Math.abs(l + r) ? l - r : 0;\n }\n return (rx, ry, qx, qy, px, py) => \n orient(px, py, rx, ry, qx, qy) || \n orient(rx, ry, qx, qy, px, py) || \n orient(qx, qy, px, py, rx, ry);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":988,"value":"md`#### Fully robust\n\nAn exact test using [robust-predicates](https://github.com/mourner/robust-predicates), a modern JavaScript port of Jonathan Shewchuk's robust predicates code.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":467,"value":"orientRobust = (await import('https://cdn.jsdelivr.net/npm/robust-predicates/+esm')).orient2d","pinned":true,"hidden":false,"mode":"ojs"},{"id":1021,"value":"md`## Acknowledgements\n\nThanks to [Jacob Rus](https://observablehq.com/@jrus) and [Philippe Rivière](https://observablehq.com/@fil) for contributions to this notebook, [Mikola Lysenko](https://twitter.com/MikolaLysenko) for fantastic lecture notes about robustness that got me fascinated with the subject, and [Jonathan Shewchuk](https://people.eecs.berkeley.edu/~jrs/) for his robust predicates code and other amazing contributions to computer science.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":417,"value":"md`## Appendix`","pinned":false,"hidden":false,"mode":"ojs"},{"id":313,"value":"fancyRange = (options) => {\n const {width, min, max, value, step, label} = Object.assign({\n width: 640,\n min: 0,\n max: 1,\n value: 0,\n step: 'any',\n label: v => html`${v.toLocaleString()}`\n }, options);\n \n const input = html``;\n const div = html`

${input}

`;\n let labelEl;\n function update() {\n div.value = +input.value; \n if (labelEl) div.removeChild(labelEl);\n labelEl = html`${label(div.value)}`;\n div.appendChild(labelEl);\n div.dispatchEvent(new CustomEvent('input'));\n }\n input.oninput = update;\n update();\n return div;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1014,"value":"import {animate} from '@jrus/misc'","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}