31:Tb9c,md`# Approximating Bézier Curves for CNC Grbl and several other CNC controllers don't support fancy [Bézier curve commands](http://linuxcnc.org/docs/2.6/html/gcode/gcode.html#sec:G5-Cubic-Spline) in their dialect of G-Code _(update Jan. 2022: I have learned that the B-spline in the linked article is *not* a Bézier curve. It's a [basis spline](https://en.wikipedia.org/wiki/B-spline).)_, but they do tend to support [linear interpolation and circular interpolation](https://wiki.shapeoko.com/index.php/G-Code#Motion_.28G.29). After learning this, I became interested in techniques for approximating Bézier curves using lines and circular arcs and wanted to compare the results. Grbl eventually interpolates circular arcs into a series of line segments, so linear interpolation _before_ sending G-Code to the controller wouldn't materially affect the precision of the toolpath assuming the curve is subdivided in small enough increments. However, using line segments results in an explosion of G1 commands to move the tool in those small lines, and this fills up Grbl's serial buffer. With no knowledge of Grbl's motion planning, I assume a filled-up buffer makes it difficult for Grbl to look ahead for future instructions and plan out its accelerations, resulting in very slow motion along curves. For a pen plotter, this means an inky pen will blot out on the paper, and that's not ideal (unless you want that kind of thing). Circular arcs can describe a curve with nice precision and generally fewer instructions to send over the wire than linear interpolation, and I wanted to see that here in this notebook. I haven't confirmed whether this improves the acceleration management on the controller yet. Grbl might be making those plans _after_ the circular arcs have been blown out into tiny line segments. I figure it's still a win if we capture an ideal geometry with fewer commands. *Quick weekend update! I learned the [difference between G93 and G94 commands](http://linuxcnc.org/docs/2.6/html/gcode/gcode.html#sec:G93-G94-G95-Mode) today. I suspect some of the behaviors I was seeing on my machine that prompted this exploration were due to my CAM software using G93 (Inverse Time Mode) where the controller attempts to complete each operation in a fixed amount of time. For operations with small distances to travel, the machine will move very slowly to fill the time. G94 (Units per Minute Mode) tells the controller to move the controlled point at a constant speed regardless of the distance per operation. Initial experiments with my plotter using G94 indicate that its behavior with linearly interpolated curves _and_ circularly approximated curves are pretty similar. I'm honestly not sure why anyone would use G93 mode given how much complexity it leaves for the end-user, so I'm curious why it exists.* *Also, G95 (Units per Revolution Mode) sounds awesome, but my CNC does not have a spindle at the moment. I'm curious how it's used though.* `32:T81a,function createBezierCurveTo_circular(tolerance, debuggin = false) { // Based on https://pomax.github.io/bezierinfo/#arcapproximation return function bezierCurveTo_circular( ctx, x1, y1, cp1x, cp1y, cp2x, cp2y, x2, y2 ) { function pointForT(t) { return bezierCurve(t, [x1, y1], [cp1x, cp1y], [cp2x, cp2y], [x2, y2]); } const colors = colorRotator(); let t = 0; while (t < 1) { let e = 1; // Yikes! Does this algorithm always stop? I think so... while (true) { const n = (t + e) / 2; const p1 = pointForT(t); const p2 = pointForT(n); const p3 = pointForT(e); const { center, radius } = circleFromThreePoints(p1, p2, p3); // Choose two points between curve segments t to n and n to e const check1 = pointForT((t + n) / 2); const check2 = pointForT((n + e) / 2); const error1 = computeCircularError(check1, center, radius); const error2 = computeCircularError(check2, center, radius); const maxError = Math.max(error1, error2); if (maxError < tolerance) { const startAngle = angle(center, p1); const middleAngle = angle(center, p2); const endAngle = angle(center, p3); ctx.beginPath(); ctx.strokeStyle = colors.next().value; ctx.arc( ...center, radius, startAngle, endAngle, // I don't understand this following logic for CW/CCW, but it works... startAngle > middleAngle || middleAngle > endAngle ); ctx.stroke(); if (debuggin) { ctx.beginPath(); ctx.globalAlpha = 0.25; ctx.lineWidth = 1; ctx.arc(...center, radius, 0, 2 * Math.PI); ctx.stroke(); ctx.lineWidth = 2; ctx.globalAlpha = 1; } t = e; break; } else { // Not quite binary search... e = n; } } } }; }33:T707,function createBezierCurveTo_circularSubdivided( subdivisions, debuggin = false ) { // Based on https://pomax.github.io/bezierinfo/#arcapproximation return function bezierCurveTo_circular( ctx, x1, y1, cp1x, cp1y, cp2x, cp2y, x2, y2 ) { function pointForT(t) { return bezierCurve(t, [x1, y1], [cp1x, cp1y], [cp2x, cp2y], [x2, y2]); } const colors = colorRotator(); for (let i = 1; i <= subdivisions; i++) { const t = (i - 1) / subdivisions; const e = i / subdivisions; const n = (t + e) / 2; const p1 = pointForT(t); const p2 = pointForT(n); const p3 = pointForT(e); const { center, radius } = circleFromThreePoints(p1, p2, p3); // Choose two points between curve segments t to n and n to e const check1 = pointForT((t + n) / 2); const check2 = pointForT((n + e) / 2); const error1 = computeCircularError(check1, center, radius); const error2 = computeCircularError(check2, center, radius); const maxError = Math.max(error1, error2); const startAngle = angle(center, p1); const middleAngle = angle(center, p2); const endAngle = angle(center, p3); ctx.beginPath(); ctx.strokeStyle = colors.next().value; ctx.arc( ...center, radius, startAngle, endAngle, // I don't understand this following logic for CW/CCW, but it works... startAngle > middleAngle || middleAngle > endAngle ); ctx.stroke(); if (debuggin) { ctx.beginPath(); ctx.globalAlpha = 0.25; ctx.lineWidth = 1; ctx.arc(...center, radius, 0, 2 * Math.PI); ctx.stroke(); ctx.lineWidth = 2; ctx.globalAlpha = 1; } } }; }34:[] 2f:{"title":"Approximating Bézier Curves for CNC","cells":[{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":370,"value":"md`# Linear Interpolation`","pinned":false,"hidden":false,"mode":"ojs"},{"id":98,"value":"makeFigure(\n width,\n width / 2,\n 'Subdividing the curve into 10 segments',\n createBezierCurveTo_subdivided(10)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":873,"value":"{\n const n = 4;\n const figWidth = width / n;\n const figHeight = width / (n + 1);\n\n return html`\n

\n ${makeFigure(\n figWidth,\n figHeight,\n '5 subdivisions',\n createBezierCurveTo_subdivided(5)\n )}\n ${makeFigure(\n figWidth,\n figHeight,\n '10 subdivisions',\n createBezierCurveTo_subdivided(10)\n )}\n ${makeFigure(\n figWidth,\n figHeight,\n '25 subdivisions',\n createBezierCurveTo_subdivided(25)\n )}\n ${makeFigure(\n figWidth,\n figHeight,\n '50 subdivisions',\n createBezierCurveTo_subdivided(50)\n )}\n

\n`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":7,"value":"function createBezierCurveTo_subdivided(subdivisions) {\n return function bezierCurveTo_subdivided(\n ctx,\n x1,\n y1,\n cp1x,\n cp1y,\n cp2x,\n cp2y,\n x2,\n y2\n ) {\n const colors = colorRotator();\n for (let i = 1; i <= subdivisions; i++) {\n const t = i / subdivisions;\n const tx = b3(t, x1, cp1x, cp2x, x2);\n const ty = b3(t, y1, cp1y, cp2y, y2);\n\n ctx.strokeStyle = colors.next().value;\n ctx.lineTo(tx, ty);\n ctx.stroke();\n ctx.beginPath();\n ctx.moveTo(tx, ty); // Wonky but easier than keeping prev coordinates...\n }\n };\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":375,"value":"md`## Using an ideal step size\n\nThis is not much different than the subdivision approach, but it results in a more consistent level of precision no matter the wiggly-ness of the curve. The number of subdivisions is determined by the ideal stepsize and the approximate length, rather than a constant number of subdivisions.\n\nThere is [a better way](https://pomax.github.io/bezierinfo/#tracing) to do this.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":184,"value":"makeFigure(\n width,\n width / 2,\n 'Using segments with approximately 100px length',\n createBezierCurveTo_stepsize(100)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":175,"value":"function createBezierCurveTo_stepsize(idealStepSize) {\n return function bezierCurveTo_stepsize(\n ctx,\n x1,\n y1,\n cp1x,\n cp1y,\n cp2x,\n cp2y,\n x2,\n y2\n ) {\n const length = arclength_approx(20, x1, y1, cp1x, cp1y, cp2x, cp2y, x2, y2);\n const subdivisions = Math.ceil(length / idealStepSize);\n const colors = colorRotator();\n\n for (let i = 1; i <= subdivisions; i++) {\n const t = i / subdivisions;\n const tx = b3(t, x1, cp1x, cp2x, x2);\n const ty = b3(t, y1, cp1y, cp2y, y2);\n\n ctx.strokeStyle = colors.next().value;\n ctx.lineTo(tx, ty);\n ctx.stroke();\n ctx.beginPath();\n ctx.moveTo(tx, ty); // Wonky but easier than keeping prev coordinates...\n }\n };\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":846,"value":"md`There's the possibility to dynamically search for an appropriately long line segment by computing the error of the line segment compared to the original curve. I haven't explored that yet, but this technique is essentially what's going on in the next section, except we're using circles rather than lines.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":395,"value":"md`# Circular Approximation\n\nStraight lines are cool, but we're trying to make bendy lines. Circles are bendy.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":479,"value":"makeFigure(width, width / 2, '', createBezierCurveTo_circular(0.25, true))","pinned":false,"hidden":false,"mode":"ojs"},{"id":620,"value":"{\n const n = 4;\n const figWidth = width / n;\n const figHeight = width / (n + 1);\n\n return html`\n

\n ${makeFigure(\n figWidth,\n figHeight,\n 'High tolerance',\n createBezierCurveTo_circular(5)\n )}\n ${makeFigure(\n figWidth,\n figHeight,\n 'Medium tolerance',\n createBezierCurveTo_circular(2.5)\n )}\n ${makeFigure(\n figWidth,\n figHeight,\n 'Low tolerance',\n createBezierCurveTo_circular(1)\n )}\n ${makeFigure(\n figWidth,\n figHeight,\n 'Tiny tolerance',\n createBezierCurveTo_circular(0.1)\n )}\n

\n`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":929,"value":"md`You can see some sharp corners on the left two examples. There's a better way to [approximate Bézier curves using biarcs](http://dlacko.org/blog/2016/10/19/approximating-bezier-curves-by-biarcs/) instead where every circular arc is guaranteed to share the same tangent at the joint as its predecessor. I didn't use that approach here, as you can see.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":405,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":415,"value":"// Compare the distance from a point on the curve to the center of the circle to the radius of the circle.\n// We want this number to be 0.\nfunction computeCircularError(pointOnCurve, center, radius) {\n const actual = distance(pointOnCurve, center);\n const expected = radius;\n return Math.abs(actual - expected);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":897,"value":"md`Just for fun, here's an implementation with a fixed number of subdivisions.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":884,"value":"makeFigure(\n width,\n width / 2,\n '10 subdivisions',\n createBezierCurveTo_circularSubdivided(10, true)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":904,"value":"{\n const n = 4;\n const figWidth = width / n;\n const figHeight = width / (n + 1);\n\n return html`\n

\n ${makeFigure(\n figWidth,\n figHeight,\n '1 subdivision',\n createBezierCurveTo_circularSubdivided(1)\n )}\n ${makeFigure(\n figWidth,\n figHeight,\n '2 subdivisions (lol)',\n createBezierCurveTo_circularSubdivided(2)\n )}\n ${makeFigure(\n figWidth,\n figHeight,\n '5 subdivisions',\n createBezierCurveTo_circularSubdivided(5)\n )}\n ${makeFigure(\n figWidth,\n figHeight,\n '10 subdivisions',\n createBezierCurveTo_circularSubdivided(10)\n )}\n

\n`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":880,"value":"$33","pinned":true,"hidden":false,"mode":"ojs"},{"id":124,"value":"md`# Helpers`","pinned":false,"hidden":false,"mode":"ojs"},{"id":703,"value":"function angle([x1, y1], [x2, y2]) {\n const dx = x2 - x1;\n const dy = y2 - y1;\n return Math.atan2(dy, dx);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":652,"value":"function makeFigure(width, height, caption, fn) {\n const offset = 30;\n\n return html`\n

\n ${compareBezier(\n width,\n height,\n [offset, offset],\n [offset, height * 2 - offset],\n [width + offset, offset],\n [width - offset, height - offset],\n fn\n )}\n

${caption}

\n \n`;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":647,"value":"function* colorRotator() {\n const colors = [\n '#FF4C6D',\n '#00D3FF',\n '#FF964E',\n '#00F1BA',\n '#DF85FF',\n '#FFCC2C'\n ];\n\n let index = 0;\n\n while (true) {\n const color = colors[index];\n index = (index + 1) % colors.length;\n yield color;\n }\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":126,"value":"// Draw a reference curve and a comparison curve using the same control points\nfunction compareBezier(width, height, start, cp1, cp2, end, pathFn) {\n const ctx = DOM.context2d(width, height, window.devicePixelRatio);\n\n // Reference curve using browser API\n ctx.beginPath();\n ctx.lineWidth = 8;\n ctx.strokeStyle = '#eee';\n ctx.moveTo(...start);\n ctx.bezierCurveTo(...cp1, ...cp2, ...end);\n ctx.stroke();\n\n // Compare curve using provided pathFn\n ctx.beginPath();\n ctx.lineWidth = 2;\n ctx.moveTo(...start);\n pathFn(ctx, ...start, ...cp1, ...cp2, ...end);\n\n return ctx.canvas;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":276,"value":"md`## Geometry`","pinned":false,"hidden":false,"mode":"ojs"},{"id":140,"value":"function distance_4(x1, y1, x2, y2) {\n return Math.sqrt(Math.pow(x2 - x1, 2) + Math.pow(y2 - y1, 2));\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":311,"value":"function distance(p1, p2) {\n return distance_4(...p1, ...p2);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":325,"value":"function midpoint([x1, y1], [x2, y2]) {\n return [(x1 + x2) / 2, (y1 + y2) / 2];\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":330,"value":"// Finds a coordinate describing a line perpendicular to the line\n// from p1 to p2 that passes through p2\nfunction perpendicularToLine([x1, y1], [x2, y2]) {\n return [\n x2 - (25 * (y1 - y2)) / distance([x1, y1], [x2, y2]),\n y2 + (25 * (x1 - x2)) / distance([x1, y1], [x2, y2])\n ];\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":327,"value":"function circleFromThreePoints(p1, p2, p3) {\n const mp1 = midpoint(p1, p2);\n const mp2 = midpoint(p1, p3);\n\n const perp1 = perpendicularToLine(p1, mp1);\n const perp2 = perpendicularToLine(p1, mp2);\n\n const center = lineLineIntersection([perp1, mp1], [perp2, mp2]);\n\n const radius = distance(center, p1);\n return { center, radius };\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":285,"value":"// From https://pomax.github.io/bezierinfo/#intersections\nfunction lineLineIntersection_8(x1, y1, x2, y2, x3, y3, x4, y4) {\n const nx = (x1 * y2 - y1 * x2) * (x3 - x4) - (x1 - x2) * (x3 * y4 - y3 * x4);\n const ny = (x1 * y2 - y1 * x2) * (y3 - y4) - (y1 - y2) * (x3 * y4 - y3 * x4);\n const d = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4);\n\n // If d === 0, the lines are paralllel and there is no intersection.\n // This will return the coordinate [Infinity, Infinity] which is suitable for use here.\n return [nx / d, ny / d];\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":282,"value":"// From https://pomax.github.io/bezierinfo/#intersections\nfunction lineLineIntersection_4(p1, p2, p3, p4) {\n const [x1, y1] = p1;\n const [x2, y2] = p2;\n const [x3, y3] = p3;\n const [x4, y4] = p4;\n return lineLineIntersection_8(x1, y1, x2, y2, x3, y3, x4, y4);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":294,"value":"// From https://pomax.github.io/bezierinfo/#intersections\nfunction lineLineIntersection(line1, line2) {\n return lineLineIntersection_4(...line1, ...line2);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":34,"value":"md`\n## Bézier Functions\n\nFrom https://en.wikipedia.org/wiki/B%C3%A9zier_curve\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":21,"value":"function b3p0(t, p) {\n var k = 1 - t;\n return k * k * k * p;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":23,"value":"function b3p1(t, p) {\n var k = 1 - t;\n return 3 * k * k * t * p;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":25,"value":"function b3p2(t, p) {\n var k = 1 - t;\n return 3 * k * t * t * p;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":27,"value":"function b3p3(t, p) {\n return t * t * t * p;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":9,"value":"// Compute value on a single dimension for bezier curve with given control points\nfunction b3(t, p0, p1, p2, p3) {\n return b3p0(t, p0) + b3p1(t, p1) + b3p2(t, p2) + b3p3(t, p3);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":421,"value":"function bezierCurve(t, [p0x, p0y], [p1x, p1y], [p2x, p2y], [p3x, p3y]) {\n return [b3(t, p0x, p1x, p2x, p3x), b3(t, p0y, p1y, p2y, p3y)];\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":384,"value":"// Approximate the length of a curve via linear interpolation. Could be very off!\nfunction arclength_approx(\n subdivisions,\n x1,\n y1,\n cp1x,\n cp1y,\n cp2x,\n cp2y,\n x2,\n y2\n) {\n let px = x1;\n let py = y1;\n let length = 0;\n for (let i = 1; i <= subdivisions; i++) {\n const t = i / subdivisions;\n const tx = b3(t, x1, cp1x, cp2x, x2);\n const ty = b3(t, y1, cp1y, cp2y, y2);\n length += distance_4(px, py, tx, ty);\n px = tx;\n py = ty;\n }\n return length;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":599,"value":"html`\n`","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}