31:T648,Explanation = md` We can use the magic of [Möbius transformations][] and the [cross-ratio][] to construct an interpolant ${$`P(t)`} tracing the circular arc from point ${$`A`} to point ${$`B`}, passing through point ${$`M`}, such that ${$`P(0) = A`}, ${$`P(\tfrac{1}{2}) = M`}, ${$`P(1) = B`}. Note that points ${$`A`}, ${$`B`}, and ${$`M`} must be distinct. [Möbius transformations]: https://en.wikipedia.org/wiki/Möbius_transformation [cross-ratio]: https://en.wikipedia.org/wiki/Cross-ratio If we want to be careful, we can separately define points and displacement vectors between points and construct “complex numbers” as the ratios of those vectors. Or we can be sloppy and pretend that planar points, planar vectors, and complex numbers are all the same type of object. We can at least give a nod to this sloppiness by labeling our complex values ${$`a`}, ${$`b`}, and ${$`m`} and our function ${$`z(t)`} instead. To construct our interpolant we equate the cross-ratios: ${$$` (z(t), a; m, b) = (t, 0; \tfrac{1}{2}, 1) `} That is: ${$$` \dfrac{(z(t) - a)(m - b)} {(z(t) - b)(m - a)} = \dfrac{(t - 0)(\tfrac{1}{2} - 1)} {(t - 1)(\tfrac{1}{2} - 0)} = \dfrac{t}{1 - t} `} In computer graphics terms we can simplify this to ${$$` z(t) = \dfrac{\operatorname{lerp}((b - m)a, (m - a)b, t)} {\operatorname{lerp}(b - m, m - a, t)} `} where all of the arguments and operators involved are complex. As desired, ${$$` \begin{aligned} z(0) &= a (b - m) / (b - m) = a \\ z(1) &= b (m - a) / (m - a) = b \\ z(\tfrac{1}{2}) &= m(b - a) / (b - a) = m \end{aligned} `} `32:T408,md` This rational quadratic is substantially cheaper than using trigonometry and both cheaper and better behaving than interpolating along a chord and then projecting onto the circle along a line through its center. It costs only 1 complex division and 2 complex [lerps][] for each evaluation. [lerps]: https://en.wikipedia.org/wiki/Linear_interpolation Importantly, this works for points in a straight line (a circle of infinite diameter) and we never need to explicitly construct the center or radius of the circle. Given points already known to lie on a particular circle, applying this method requires first computing ${$`m`}, the midpoint of the great-circle arc from ${$`a`} to ${$`b`}, which costs one square root and one division. For every subsequent evaluation we need only one division. For moderate angles with the midpoint computed precisely, this parametrization gets very close to arclength-uniform; even for a half-circle arc, the speed at the middle of the arc is only twice as fast as the speed at the ends. `33:Tb29,Noncommutative = md` ## Generalization for noncommutative values This method can also be useful for interpolating on a circular arc in ${$`n`}-dimensional space, with values which might be Euclidean points or (multi-)vectors: for example points on a 2-sphere embedded in 3-dimensional space or unit quaternions (which lie on a 3-sphere). But in cases where the values may not be commutative, we need to be a bit more careful than we were above. Remember, for commutative values we can write (letting ${$`\tilde{t} = 1 - t`}): ${$$` z(t) = \dfrac{(b - m)a\tilde{t} + (m - a)bt} {(b - m)\tilde{t} + (m - a)t} `} For noncommutative values we instead need to put the denominator ${$`q(t) = (b - m)\tilde{t} + (m - a)t`} into the middle of each term separately: ${$$` x(t) = (b - m)q(t)^{-1} a \tilde{t} + (m - a)q(t)^{-1} b t `} We still have, as desired: ${$$` \begin{aligned} x(0) &= (b - m)(b - m)^{-1}a = a \\[.3em] x(1) &= (m - a)(m - a)^{-1}b = b \\[.3em] x\big(\tfrac{1}{2}\big) &= (b - m)(b - a)^{-1}a + (m - a)(b - a)^{-1}b \\[.1em] &= m(b - a)^{-1}(b - a) = m \end{aligned}`} The last of these uses the identity ${$`b(b - a)^{-1}a = \left(a^{-1} - b^{-1}\right)^{-1} = a(b - a)^{-1}b`}. The above, as written, takes a bit more computation than necessary, as it requires more multivector multiplication than before, but is still relatively quite a cheap operation. If we prefer we can write this asymmetrically as: ${$$` \begin{aligned} x(t) &= a + (m - a)q(t)^{-1} (b - a) t \\[.1em] &= b - (b - m)q(t)^{-1} (b - a) \tilde{t} \end{aligned}`} For practical computation, if we assume ${$`a`}, ${$`m`}, and ${$`b`} are Euclidean vectors, and multiplication is the [geometric product](https://en.wikipedia.org/wiki/Geometric_algebra), we can first define ${$`u = m - a`}, ${$`v = b - m`}, ${$`h = b - a`}, and then express our interpolating function as: ${$$` x(t) = a + \dfrac{t}{\operatorname{lerp}(v, u, t)^2}\operatorname{lerp}(uvh, u^2h, t) `} Or if the values are more general invertible multivectors, we can use the multivector reverse ${$`w^\dagger`}, which generalizes the complex conjugate, ${$`w^{-1} = w^\dagger / (w w^\dagger) = w^\dagger / |w|^2`}: ${$$` \begin{aligned} x(t) &= a + \dfrac{t}{|\operatorname{lerp}(v, u, t)|^2}\operatorname{lerp}(uv^\dagger h, |u|^2h, t) \\[1em] &= b + \dfrac{t - 1}{|\operatorname{lerp}(v, u, t)|^2}\operatorname{lerp}(|v|^2h, vu^\dagger h, t) \end{aligned}`} The initial setup cost is 3 vector subtractions, 2 multivector–multivector multiplications, 1 dot product, and 1 scalar–vector multiplication. The per-point runtime cost is 2 vector lerps, 1 dot product, 1 scalar division, 1 scalar–vector multiplication, and 1 vector addition. ` 34:T5dc,plotarc = () => { let update; const height = 500, p1 = [width / 5, 200], p2 = [width / 3, 100], p3 = [2 * width / 3, 350], svg = d3.select(DOM.svg(width, height)), arc = svg.append("path") .attr("fill", "none") .attr("stroke-width", "2px") .attr("stroke", "#79e"), small_points = svg.append("g") .selectAll("circle") .data(Array.from(Array(33), (_,i) => i/32)) .enter().append("circle") .attr("fill", "#c38") .attr("r", 3), point = svg.append("g") .attr("cursor", "move") .attr("pointer-events", "all") .attr("stroke", "transparent") .attr("stroke-width", 30) .selectAll("circle") .data([p1, p2, p3]) .enter().append("circle") .attr("r", 7) .attr("fill", "#26b") .call(d3.drag() .subject(([x, y]) => ({x, y})) .on("drag", (d) => { d[0] = clamp(5, width-5, d3.event.x); d[1] = clamp(5, height-5, d3.event.y); update(); })); update = function update() { const parametric_arc = circle_interp(p1, p2, p3), arc_path = bezpts_to_svgpath(bezeval( parametric_arc, [0,1], [0,width, 0, height], 0.1)); arc.attr("d", arc_path); point .attr("cx", d => d[0]) .attr("cy", d => d[1]); small_points .attr("cx", t => parametric_arc(t)[0]) .attr("cy", t => parametric_arc(t)[1]); } update(); return svg.node(); }35:[] 2f:{"title":"(Rational) Circular Arc Interpolation","cells":[{"id":0,"value":"Diagram = md`# (Rational) Circular Arc Interpolation`","pinned":false,"hidden":false,"mode":"ojs"},{"id":143,"value":"plotarc()","pinned":false,"hidden":false,"mode":"ojs"},{"id":5,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":125,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":207,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":333,"value":"Infty = md`\n## Dealing with infinities\n\nIf any of ${$`a`}, ${$`b`}, or ${$`m`} is infinite, our formulas above turn into ${$`\\infty/\\infty`}, which is not very practically useful. But we can factor the common ${$`\\infty`} out of numerator and denominator to arrive at simplified special cases.\n\nIf ${$`b = \\infty`}:\n\n${$$`\nz(t) = a + \\dfrac{t}{1 - t}(m - a)\n`}\n\nIf ${$`a = \\infty`}:\n\n${$$`\nz(t) = b + \\dfrac{1 - t}{t}(m - b)\n`}\n\n\nFinally, if ${$`m = \\infty`}:\n\n${$$`\nz(t) = a + \\dfrac{t}{2t - 1}(b - a)\n`}\n\nIf ${$`a`}, ${$`b`}, and ${$`m`} are co-linear, and ${$`m`} is not between ${$`a`} and ${$`b`}, then ${$`z(t)`} will be infinite for some value of ${$`t \\in [0, 1]`}. This is expected (and often desirable) behavior.\n\n* * *\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":58,"value":"import { csub, cmul, cdiv } from \"39dbf7de6ad6c799\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":56,"value":"clerp = ([x0, y0], [x1, y1], t) =>\n [x0*(1-t) + x1*t, y0*(1-t) + y1*t]","pinned":true,"hidden":false,"mode":"ojs"},{"id":55,"value":"circle_interp = (a, m, b) => {\n const\n b_m = csub(b, m), \n m_a = csub(m,a),\n ab_m = cmul(a, b_m),\n bm_a = cmul(b, m_a);\n return (t) => cdiv(\n clerp(ab_m, bm_a, t),\n clerp(b_m, m_a, t));\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":29,"value":"$34","pinned":true,"hidden":false,"mode":"ojs"},{"id":198,"value":"clamp = (a, b, x) => Math.max(a, Math.min(x, b))","pinned":true,"hidden":false,"mode":"ojs"},{"id":66,"value":"import {bezeval, bezpts_to_svgpath} from '@jrus/bezplot'","pinned":false,"hidden":false,"mode":"ojs"},{"id":6,"value":"import {math_css, $, $$} with {$css as $css} from \"@jrus/misc\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":9,"value":"$$css = html`${math_css}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":32,"value":"d3 = require('d3@5')","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1","titling":"auto"}