31:T774,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\bigl(\tfrac{1}{2}\bigr) = 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) = \left( t, 0; \tfrac{1}{2}, 1 \right) `} That is: ${$$` \dfrac{(z(t) - a)(m - b)} {(z(t) - b)(m - a)} = \dfrac{(t - 0)\bigl( \tfrac{1}{2} - 1 \bigr)} {(t - 1)\bigl( \tfrac{1}{2} - 0 \bigr)} = \dfrac{t}{1 - t} `} In computer graphics terms we can simplify this to ${$$`\boxed{\; z(t) = \dfrac{\operatorname{lerp}((b - m)a, (m - a)b, t)} {\operatorname{lerp}(b - m, m - a, t)} \;} \tag{1} `} where all of the arguments and operators involved are complex. As desired, ${$$` \begin{aligned} z(0) &= (b - m) a \big/ (b - m) = a \\[.3em] z(\tfrac{1}{2}) &= (b - a) m \big/ (b - a) = m \\[.3em] z(1) &= (m - a) b \big/ (m - a) = b \end{aligned} `} If ${$`m`} happens to be the straight-line midpoint ${$`m = \tfrac{1}{2}(a + b)`}, then ${$`m - a = b - m`} and this formula degenerates to linear interpolation, ${$`z(t) = \operatorname{lerp}(a, b, t)`}. `32:T571,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 We can somewhat see why this is so nice by looking at the Taylor expansion for the half-angle tangent. For relatively small angles, the parametrization is nearly arclength-linear: ${$$`\begin{aligned} 2\tan\tfrac12 \theta &= \theta + \tfrac1{12}\theta^3 + \mathcal{O}(\theta^5) \\[10mu] 2\arctan\tfrac12 x &= x - \tfrac1{12}x^3 + \mathcal{O}(x^5) \end{aligned}`} 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 arc from ${$`a`} to ${$`b`}, which typically 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:T4be,Infty = md` ## Dealing with infinities If 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. Notice that ${$`\operatorname{lerp}(0, 1, t) = t`}, ${$`\operatorname{lerp}(1, 0, t) = 1 - t`}, and ${$`\operatorname{lerp}(-1, 1, t) = 2t - 1`}: If ${$`b = \infty`}: ${$$`\begin{aligned} z(t) &= a + \dfrac{t}{1 - t}(m - a) \\[1.5em] z'(t) &= \dfrac{m - a}{(1 - t)^2} \end{aligned}`} If ${$`a = \infty`}: ${$$`\begin{aligned} z(t) &= b - \dfrac{1 - t}{t}(b - m) \\[1.5em] z'(t) &= \dfrac{b - m}{t^2} \end{aligned}`} Finally, if ${$`m = \infty`}: ${$$`\begin{aligned} z(t) &= a + \dfrac{t}{2t - 1}(b - a) \\[1.5em] z'(t) &= - \dfrac{b - a}{(2t - 1)^2} \end{aligned}`} ${$$` `} If ${$`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. `34:Tb31,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{1 - t}{|\operatorname{lerp}(v, u, t)|^2}\operatorname{lerp}(|v|^2h, vu^\dagger h, t) \tag{5} \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. ` 35: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(); }36:[] 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":377,"value":"Derivative = md`## Derivative\n\nWe can use the quotient rule to find the derivative of (1), which simplifies nicely:\n\n${$$`\nz'(t) = \\dfrac{(b - a)(b - m)(m - a)}\n {\\operatorname{lerp}(b - m, m - a, t)^2}\n\\tag{2}\n`}\n\nAt the key points,\n\n${$$` \\begin{aligned}\nz'(0) &= (b - a)(m - a) / (b - m) \\\\[.3em]\nz'(\\tfrac{1}{2}) &= 4(b - m)(m - a) / (b - a) \\\\[.3em]\nz'(1) &= (b - a)(b - m) / (m - a)\n\\end{aligned}\n`}\n\n\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":398,"value":"md`\n## Starting from the origin\n\nWhen ${$`a = 0`}, our formula gets even simpler:\n\n${$$`\\begin{aligned}\nz(t) &= \\dfrac{mbt}{\\operatorname{lerp}(b - m, m, t)} \\\\[1.5em]\nz'(t) &= \\dfrac{mb(b - m)}{\\operatorname{lerp}(b - m, m, t)^2}\n\\end{aligned}`}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":333,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":425,"value":"md`\n## Inverse\n\nSometimes it is helpful, when given a point on the arc, to figure out its parameter value (for a point off the arc, the parameter will be complex). If we let ${$`r = \\left(m - a\\right) / \\left(b - m\\right)`}, we have:\n\n${$$`\nt(z) = \\dfrac{z - a}\n {(z - a) - r(z - b)}\n\\tag{3}\n`}\n\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":207,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":470,"value":"md`\n## Prior work\n\nThe basic content of this notebook is described in Severinas Zubė (2006) [“A circle representation using complex and quaternion numbers”](https://doi.org/10.1007/s10986-006-0025-z), *Lithuanian Mathematical Journal* 46, pp. 246–255. Except Zubė uses arbitrary complex numbers as weights instead of ${$`m-a`} and ${$`b-m`}, and therefore does not notice the point ${`m`} lying on the arc.\n\nThere is also a similar idea expressed in Xuexiang Li & Junxiao Xue (2014) [“Circular Arc Representation based on Parametric Complex Functions”](https://doi.org/10.2991/lemcs-14.2014.107), LEMCS 2014\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":393,"value":"md`* * *`","pinned":false,"hidden":false,"mode":"ojs"},{"id":58,"value":"import { csub, cmul, cdiv } from '@jrus/complex'","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":"$35","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":"$W36","stdlib":"1","titling":"auto"}