{ "metadata": { "language": "Julia", "name": "", "signature": "sha256:a8f4868cde474aef4d80e3967c3f7bfab45ab2d8570ed98eb4fc75d20f4aced6" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Interpolating a function via circular arcs in Julia " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Time for some recreational math.\n", "\n", "Let's try and write a program, which given a function equation and a segment, for example:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$ f(x) = 6x^2 + 4x - 8, x \\in [-4, 4] $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "would build an approximation of this with *a sequence of circular arcs*.\n", "\n", "As a bonus, we would like to generate the corresponding [G-codes](http://en.wikipedia.org/wiki/G-code) for these arcs, so the output of our program could be fed to some automatic cutting tool, which would cut out a nice parabola, so we could make a pendant, put it on the neck and proudly wear as a proof of our accomplishments!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Possible solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Given that the problem has real-world practical applications (in 3D printing, machining et al), there has been a wealth of research on the topic.\n", "\n", "I've taken [this ancient thesis](https://smartech.gatech.edu/handle/1853/17818) as a starting point, mostly because of the simplicity of the approach it is taking.\n", "\n", "So what we are going to try is essentially a rather brute force, [greedy](http://en.wikipedia.org/wiki/Greedy_algorithm) algorithm.\n", "\n", "For some given segment **[x1, x2]**, we draw the arc approximation as a circle, that passes through the point **(x1, f(x1))**, the point **(x2, f(x2))** and also the third, middle point **(xm, f(xm))**, which for simplicity we can assume is equidistant from x1 and x2:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![](circle_3_points.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then the algorithm proceeds as following:\n", "* Start from trying to fit an arc to the whole segment **[a, b]** (so that **x1**=**a** and **x2**=**b**) and compute the error value\n", "* If error is good enough then we are done, otherwise we shrink the current segment by some ratio (say, 0.8), such that point **x1** stays in place and point **x2** moves left, repeating this until the error is good enough\n", "* Once the current arc is fit, we repeat the same thing with the next one, this time trying to cover the remaining part of the segment (i.e. **x1**=**x2_prev** and **x2**=**b**)\n", "* The previous two steps get repeated until we cover the whole segment (so that **x2**=**b** and the error is good enough)\n", "\n", "Essentially, we will try packing the arcs left-to-right, in a \"greedy\" way, until the whole segment is covered." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Julia language" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Why using Julia programming language in particular?.. [This](http://julialang.org/blog/2012/02/why-we-created-julia/) sums it up nicely.\n", "\n", "Matlab has established itself as a de-facto standard tool for scientific computations. I've been using it almost daily at the job, and it's rather good at what it does. \n", "\n", "It also has flaws:\n", "* Too expensive\n", "* Too slow, in general \n", "* The programming language itself is not very expressive\n", "\n", "The first problem is what [Octave](https://www.gnu.org/software/octave/) is trying to solve, by being free and trying to replicate most of the Matlab's functionality. However, it's even \"slower\" than Matlab, and the language itself is pretty much the same, except of a few improvements (which *was* their goal, to start with).\n", "\n", "The second problem can be somewhat mitigated in Matlab itself via vectorization (i.e. operate on the level of whole vectors/matrices as opposed to separate elements) or re-implementing performance-critical parts in C/C++.\n", "\n", "The third one is something that just had to be lived with - Matlab was apparently designed as BASIC for engineers, as simple and straightforward imperative/array programming language as possible. Which is both a benefit and a flaw. \n", "\n", "Julia is trying to solve the mentioned problems, and some others. At the moment of writing the current stable version is 0.35 (and 0.4 is being in the works), which suggests that there is still a way to go until it can be considered a mature and stable tool.\n", "\n", "But even now it seems to be doing quite a good job.\n", "\n", "And unlike Matlab, it's fun to program with, especially using a pimped REPL like [IJulia](https://github.com/JuliaLang/IJulia.jl). In fact, this article is just an IJulia notebook, which is available [here](http://nbviewer.ipython.org/github/silverio/GCode/blob/master/GCode.ipynb). So let's go step-by-step and implement the building blocks that we need for the algorithm to run." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Circle through three points" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The basic operation we need first is finding a circle that goes through the three given points.\n", "There does noes not seem a corresponding library/function available in Julia yet, so let's roll our own, having looked into the [Wikipedia article](http://en.wikipedia.org/wiki/Circumscribed_circle)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![](circumscribed_circle.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Given any three points **P1**, **P2** and **P3** we can find the circle with center **Pc** and radius **r**, circumscribed over these points using the following formulas (this is just an one of the possible forms of equations): " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$ P_c=-\\frac{\\alpha P_1 + \\beta P_2 + \\gamma P_3}{2d^2} $$\n", "$$ r=\\frac{|P_1-P_2||P_2-P_3||P_3-P_1|}{2d} $$\n", "$$d=|(P_1-P_2)\\times(P_2-P_3)|$$\n", "$$\\alpha=|P_2-P_3|^2(P_1-P_2)\\cdot(P_3-P_1)$$\n", "$$\\beta=|P_3-P_1|^2(P_1-P_2)\\cdot(P_2-P_3)$$\n", "$$\\gamma=|P_1-P_2|^2(P_3-P_1)\\cdot(P_2-P_3)$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Unit testing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But before writing any Julia code let's first figure out how could we test it.\n", "\n", "There is a Base.Test Julia module which we can include and perform some basic unit testing:" ] }, { "cell_type": "code", "collapsed": true, "input": [ "using Base.Test\n", "\n", "@test 1 + 2 == 3" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The test passes fine, there is no complaints. Let's try another one:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "@test .1 + .2 == .3" ], "language": "python", "metadata": {}, "outputs": [ { "ename": "LoadError", "evalue": "test failed: 0.1 + 0.2 == 0.3\nwhile loading In[2], in expression starting on line 1", "output_type": "pyerr", "traceback": [ "test failed: 0.1 + 0.2 == 0.3\nwhile loading In[2], in expression starting on line 1", "", " in error at error.jl:21 (repeats 2 times)" ] } ], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, this one fails. Why?.. Well, simply because [0.1 plus 0.2 is not equal 0.3](http://blog.ruslans.com/2014/12/a-story-about-angry-carrot-and-floating.html).\n", "\n", "Here's a better way:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "@test_approx_eq (.1 + .2) .3" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now it passes alright. What if we wanted to compare [tuples](http://en.wikipedia.org/wiki/Tuple)?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "@test_approx_eq (1 + 2, .1 + .2) (3, .3)" ], "language": "python", "metadata": {}, "outputs": [ { "ename": "LoadError", "evalue": "`full` has no method matching full(::(Int64,Float64))\nwhile loading In[4], in expression starting on line 1", "output_type": "pyerr", "traceback": [ "`full` has no method matching full(::(Int64,Float64))\nwhile loading In[4], in expression starting on line 1", "", " in test_approx_eq at test.jl:88", " in test_approx_eq at test.jl:119" ] } ], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": {}, "source": [ "This does not work, unfortunately, but Julia's flexible enough, so we can fix it by rolling our own unit testing macro, with Blackjack and stuff:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "is_close(x::Number, y::Number, epsilon = eps()) = \n", " abs(y - x) <= max(1, abs(x), abs(y))*epsilon\n", "\n", "is_close(x, y, epsilon = eps()) = \n", " all(x -> is_close(x..., epsilon), zip(x, y)) \n", "\n", "function test_close(a, b, eps, astr, bstr)\n", " if !is_close(a, b, eps)\n", " println(\"test failed: $(astr) is not close to $(bstr);\\n $(a) vs $(b)\")\n", " end\n", "end\n", "\n", "macro test_close(a, b)\n", " :( test_close($(a), $(b), 1e-6, $(string(a)), $(string(b))) )\n", "end\n", "\n", "@test_close (1 + 2, .1 + .2) (3, .3)\n", "@test_close (1 + 2, .1 + .2) (2, .3)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "test failed: (1 + 2,0.1 + 0.2) is not close to (2,0.3);\n", " (3,0.30000000000000004) vs (2,0.3)" ] }, { "output_type": "stream", "stream": "stdout", "text": [ "\n" ] } ], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is a bunch of language concepts already packed in a few lines here, such as function overloading/specialization, default parameter values, type inference, iterators, lambdas, tuple splicing, string interpolation, macros and what have you.\n", "\n", "But this is just a glimpse into the Julia's capabilities, so if it does not make sense to you, then let's just pretend that Julia's test library already has it and move on." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Circle through three points - the code" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Armed with this, we can finally write some code, related to the problem at hand (and the corresponding unit tests). It is a rather straightforward translation of the formulas:" ] }, { "cell_type": "code", "collapsed": true, "input": [ "# returns the (center, radius) of the circle, circumscribed around three points\n", "function circumscribed_circle(P1, P2, P3)\n", " d12 = vec(P1 - P2); n12 = dot(d12, d12)\n", " d23 = vec(P2 - P3); n23 = dot(d23, d23)\n", " d31 = vec(P3 - P1); n31 = dot(d31, d31)\n", " d = norm(cross([d12; 0], [d23; 0]))\n", " r = sqrt(n12*n23*n31)/(2d)\n", " alpha = n23*dot(d12, d31)\n", " beta = n31*dot(d12, d23)\n", " gamma = n12*dot(d31, d23)\n", " Pc = -(alpha*P1 + beta*P2 + gamma*P3)/(2d^2)\n", " (Pc, r)\n", "end\n", "\n", "# the circumscribed_circle tests\n", "@test_close ([-1, 1], 1) circumscribed_circle([-2; 1], [-1; 2], [0; 1])\n", "@test_close ([0, 0], 2) circumscribed_circle([2; 0], [0; 2], [sqrt(2); sqrt(2)])" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 6 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Point orientation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another function we'll need is the one that finds out whether three points are winded counter-clockwise (CW) or clock-wise.\n", "\n", "We can do that by looking at the sign of the third component of a vector cross product between two vectors spanning our three points. \n", "\n", "As a bonus we can also determine if all the three points are lined up (i.e they are neither CW nor CCW oriented)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![](point_orientation.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$c_z=((P2-P1)\\times(P3-P1)) \\cdot [0,0,1]$$" ] }, { "cell_type": "code", "collapsed": true, "input": [ "const EPSILON = 1e-5\n", "\n", "# finds the orientation of three 2D points\n", "# -1 for CW, 1 for CCW, 0 if on the same line \n", "function orientation(P1, P2, P3) \n", " _, _, cz = cross([vec(P2 - P1); 0], [vec(P3 - P1); 0]) \n", " if abs(cz) < EPSILON\n", " cz = 0\n", " end\n", " sign(cz)\n", "end\n", "\n", "@test 1 == orientation([0; 0], [1; 0], [2; 2])\n", "@test 0 == orientation([0; 0], [1; 1], [2; 2])\n", "@test -1 == orientation([0; 0], [2; 2], [1; 0])" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Handling three points on the same line" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we've got our circumscribed circle function, let's try some other point combination:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "circumscribed_circle([1; 1], [2; 2], [3; 3])" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 8, "text": [ "([NaN,NaN],Inf)" ] } ], "prompt_number": 8 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Oops.\n", "\n", "Apparently, the part with \"any\" three points was not quite correct: we can't make an arc through them if they happen to be on the same line. Since the problem constraint was to use only arc primitives (even though in practice G-codes do allow straight line segments as well), we'd need to remedy it somehow.\n", "\n", "A solution is to offset the middle point somewhat in the direction, orthogonal to the line in question:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![](point_offset.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The corresponding formulas and the code:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$P_{offs}=P+\\epsilon \\frac{D_{\\bot}}{|D_{\\bot}|}, D_{\\bot}=\\begin{bmatrix}-D_y\\\\ D_x\\\\\\end{bmatrix}, D = P_2 - P$$" ] }, { "cell_type": "code", "collapsed": true, "input": [ "function offset_orthogonal(P2, P3)\n", " D = P2 - P3\n", " Dp = [-D[2]; D[1]]\n", " P2 + EPSILON*Dp/norm(Dp)\n", "end\n", "\n", "@test 0 != orientation([0; 0], offset_orthogonal([1; 1], [2; 2]), [2; 2])" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 9 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Arc equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What we are going to do, essentially, is to split the function's range into segments and on each segment approximate it with another function (of an arc). The arc function looks like following (coming from the circle equation): " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$ f_{arc}(x)=c_y \\pm \\sqrt{R^2 - (x - c_x)^2}, x \\in [c_x - R, c_x + R] $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are two possibilities for the arc function: it can be either convex or concave, which corresponds to the +/- in the equation:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![](arc_equation.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The code:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "function f_arc(x, c, R, convex = true) \n", " cx, cy = c\n", " plusminus = convex ? 1 : -1\n", " cy + plusminus*sqrt(max(0, R^2 - (x - cx)^2))\n", "end\n", "\n", "@test_close 3 f_arc(1, [1; 1], 2, true)\n", "@test_close -1 f_arc(1, [1; 1], 2, false)\n", "\n", "@test_close 1 f_arc(3, [1; 1], 2, true)\n", "@test_close 1 f_arc(3, [1; 1], 2, false)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 10 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Arc data structure" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To describe an arc segment it's convenient to group several attributes into a dedicated data structure. In Julia it's defined as following:" ] }, { "cell_type": "code", "collapsed": true, "input": [ "type ArcSeg\n", " c::Array{Float64, 1} # arc center\n", " R::Float64 # arc radius\n", " a::Float64 # segment start\n", " b::Float64 # segment end\n", " convex::Bool # true if top part of the circle (otherwise concave)\n", " m::Float64 # middle point ratio\n", "end" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 11 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have field data types specified here explicitly, but generally Julia does not require that." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Making an arc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's a function to actually draw an arc through two corner points and the middle one, given the function and the segment:" ] }, { "cell_type": "code", "collapsed": true, "input": [ "# makes an arc going through the segment ends\n", "function make_arc(f, seg, middle_ratio = 0.5)\n", " a, b = seg \n", " m = a + (b - a)*middle_ratio\n", " pa, pb, pm = [a; f(a)], [b; f(b)], [m; f(m)]\n", " # handle the case when all three points are on the same line\n", " porient = orientation(pa, pm, pb)\n", " if porient == 0\n", " # points are on the same line\n", " pm = offset_orthogonal(pm, pb)\n", " end\n", " c, R = circumscribed_circle(pa, pm, pb)\n", " ArcSeg(c, R, a, b, porient < 0, middle_ratio)\n", "end\n", "\n", "# inscribe an arc into itself\n", "_arc = make_arc(x -> f_arc(x, [1; 1], 5, true), (3, 4), 0.2)\n", "@test_close [1; 1] _arc.c\n", "@test_close (5, 3, 4, 0.2) (_arc.R, _arc.a, _arc.b, _arc.m)\n", "@test _arc.convex" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 12 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Error metric" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For measuring how far away is our arc from the actual function, we'll use the so-called L1-metric distance between functions:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$E(x_1, x_2)=\\int_{x_1}^{x_2} |f(x)-f_{arc}(x)| dx$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Essentially, it's the area of undesired \"slack\" between two curves:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![](error_metric.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To find the total error on the segment, we need to integrate this distance function (which looks like camel humps in our case):" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![](error_metric_function.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finding this integral can be done in different ways, but we'll use a rather simple one, [trapezoidal integration](http://en.wikipedia.org/wiki/Trapezoidal_rule), whereas the segment is split into a number of chunks of equal width, and the total area is computed as the sum of the corresponding trapezoids:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![](trapezoidal_integration.png)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "const NUM_INTEGRATION_SAMPLES = 1000\n", "function integrate_trapezoidal(f, x1, x2, N = NUM_INTEGRATION_SAMPLES) \n", " @assert(N >= 1)\n", " dx = (x2 - x1)/N;\n", " s0 = (f(x1) + f(x2))/2\n", " (s0 + sum([f(clamp(x1 + i*dx, x1, x2)) for i in 1:(N - 1)]))*(x2 - x1)/N\n", "end\n", "\n", "# integrate parabola \n", "@test_close 2^3/3 integrate_trapezoidal(x -> x^2, 0, 2)\n", "@test_close 2^3/3 integrate_trapezoidal(x -> x^2, 0, 2, 10000)\n", "@test_close 2^3/3 integrate_trapezoidal(x -> x^2, 0, 2, 1000000)\n", "\n", "# integrate half a circle (offset by 1 upwards from X axis)\n", "@test_close (pi*3^2/2 + 3*2) integrate_trapezoidal(x -> f_arc(x, [1; 1], 3, true), 1 - 3, 1 + 3, 10000)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 13 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then the error function itself is:" ] }, { "cell_type": "code", "collapsed": true, "input": [ "function err(f, arc::ArcSeg)\n", " err_metric(x) = abs(f(x) - f_arc(x, arc.c, arc.R, arc.convex))\n", " res = integrate_trapezoidal(err_metric, arc.a, arc.b)\n", "end\n", "\n", "f(x) = 2\n", "@test is_close(2 - pi/2, err(f, ArcSeg([1; 1], 1, 0, 2, true, 0.1)), 1e-4)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 14 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that package \"Calculus\" does have some numerical integration routines, which we could theoretically use, and that's what I tried at first.\n", "\n", "That did not work out very well, though. The default integration method (Simpson's) does not seem to converge for some functions.\n", "\n", "Another method there is Monte-Carlo integration:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "using Calculus\n", "function err_old(f, arc::ArcSeg)\n", " err_metric(x) = abs(f(x) - f_arc(x, arc.c, arc.R, arc.convex))\n", " res = Calculus.monte_carlo(err_metric, arc.a, arc.b, NUM_INTEGRATION_SAMPLES)\n", "end\n", "\n", "@test_close (2 - pi/2) err_old(f, ArcSeg([1; 1], 1, 0, 2, true, 0.1))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "test failed: 2 - pi / 2 is not close to err_old(f,ArcSeg([1,1],1,0,2,true,0.1));\n", " 0.42920367320510344 vs 0.4484735688613066" ] }, { "output_type": "stream", "stream": "stdout", "text": [ "\n" ] } ], "prompt_number": 15 }, { "cell_type": "code", "collapsed": false, "input": [ "@test_close (2 - pi/2) err_old(f, ArcSeg([1; 1], 1, 0, 2, true, 0.1))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "test failed: 2 - pi / 2 is not close to err_old(f,ArcSeg([1,1],1,0,2,true,0.1));\n", " 0.42920367320510344 vs 0.4610385490322167\n" ] } ], "prompt_number": 16 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The same test *both* fails and gives different answers in different runs.\n", "\n", "Non-deterministic == not good, so let's stick to our own trapezoidal wheel we've just reinvented." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The \"ears\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, our error metric is still not good enough, since it does not account for certain corner cases (\"ears\") such as here:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![](error_metric_ears1.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The orange area is the one that we've just computed, but the pink and the green \"ears\" are the extra pieces that also add to the error value.\n", "\n", "Since we generally take three arbitrary points to draw a circle through, such point configuration is not uncommon, and we need to properly adjust the error metric to account for it. \n", "\n", "Wikipedia is, again, to the resque for the [ears area](http://en.wikipedia.org/wiki/Circular_segment):" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$S=\\frac{R^2}{2}(\\theta - \\sin{\\theta}), \\theta=2\\arccos(1-\\frac{x}{R})$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![](circle_segment_area.png)" ] }, { "cell_type": "code", "collapsed": true, "input": [ "function circle_segment_area(x, R)\n", " theta = 2acos((R-x)/R)\n", " R^2*(theta - sin(theta))/2\n", "end\n", "\n", "@test_approx_eq circle_segment_area(0, 7) 0\n", "@test_approx_eq circle_segment_area(14, 7) pi*7^2\n", "@test_approx_eq circle_segment_area(7, 7) pi*7^2/2\n", "@test_approx_eq circle_segment_area(2, 7) (pi*7^2 - circle_segment_area(7*2 - 2, 7))" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 17 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Full error function" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Having taken the ears into account, let's write the full error function:" ] }, { "cell_type": "code", "collapsed": true, "input": [ "function err1(f, arc::ArcSeg)\n", " res = err(f, arc)\n", " # add the \"ears\" error\n", " fa, fb = f(arc.a), f(arc.b)\n", " cx, cy = arc.c[1], arc.c[2]\n", " if (fa > cy) != arc.convex\n", " # left \"ear\"\n", " res += circle_segment_area(arc.a - (cx - arc.R), arc.R)\n", " end\n", " if (fb > cy) != arc.convex\n", " # right \"ear\"\n", " res += circle_segment_area(cx + arc.R - arc.b, arc.R)\n", " end\n", " res\n", "end\n", "\n", "f(x) = 2\n", "@test is_close(2 - pi/2, err1(f, ArcSeg([1; 1], 1, 0, 2, true, 0.1)), 1e-4)\n", "\n", "arc = ArcSeg([1; 1], 1, 0.5, 1.5, true, 0.1)\n", "d = f_arc(0.5, arc.c, arc.R)\n", "f(x) = d\n", "@test is_close(pi - circle_segment_area(d, 1), err1(f, arc), 1e-4)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 18 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fitting a sequence of arcs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we put together all the building blocks into the arc fitting procedure, following the steps discussed earlier:" ] }, { "cell_type": "code", "collapsed": true, "input": [ "const STEP_RATIO = 0.8 \n", "const pick_center_pt = (f, a, b) -> 0.5\n", "\n", "function fit_arcs(f, range, tolerance, pick_middle_pt = pick_center_pt)\n", " ra, rb = range\n", " D = rb - ra\n", " rel_tol = tolerance/D\n", " ca = ra\n", " arcs = []\n", " while ca < rb\n", " cb = rb\n", " while true\n", " middle_ratio = pick_middle_pt(f, ca, cb)\n", " arc = make_arc(f, (ca, cb), middle_ratio)\n", " E = err1(f, arc)\n", " if E <= rel_tol*(cb - ca)\n", " arcs = [arcs, arc]\n", " ca = cb\n", " break\n", " end\n", " cb = ca + STEP_RATIO*(cb - ca)\n", " end\n", " end\n", " arcs\n", "end" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 19, "text": [ "fit_arcs (generic function with 2 methods)" ] } ], "prompt_number": 19 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We've left \"pick the middle point\" as a callback, which by default simply takes the point halfway between the segment ends." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Inspecting the result" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We did not write a unit test for this one, since it's a bit too complex of a unit already, but let's try and see how we could test it otherwise.\n", "\n", "Let's take a simple function as an example:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$f(x)=\\frac{x^3}{27} + \\frac{x}{5} - 1, x \\in [-5, 4] $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![](test_func.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Julia's DataFrames library can be used to organize/display the data in tabular form (and do much more than that):" ] }, { "cell_type": "code", "collapsed": true, "input": [ "using DataFrames\n", "\n", "fmt(x::Float64) = round(x, 3)\n", "fmt(x::Array{Float64,1}) = map(fmt, x)\n", "fmt(x::Bool) = x ? \"yes\" : \"no\"\n", "fmt(x) = x\n", "\n", "function to_table(arcs)\n", " DataFrame(\n", " c = map(x->fmt(x.c), arcs), \n", " R = map(x->fmt(x.R), arcs),\n", " a = map(x->fmt(x.a), arcs), \n", " b = map(x->fmt(x.b), arcs), \n", " m = map(x->fmt(x.m), arcs), \n", " convex = map(x->fmt(x.convex), arcs))\n", "end\n", "\n", "to_table(fit_arcs(x -> x^3/27 + x/5 - 1, [-5 4], 0.1))" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
c | R | a | b | m | convex | |
---|---|---|---|---|---|---|
1 | [11.278,-11.843] | 17.092 | -5.0 | -3.792 | 0.5 | yes |
2 | [1.386,-6.401] | 5.805 | -3.792 | -2.158 | 0.5 | yes |
3 | [0.554,-5.846] | 4.868 | -2.158 | -0.14 | 0.5 | yes |
4 | [-0.888,4.714] | 5.79 | -0.14 | 1.98 | 0.5 | no |
5 | [-0.935,3.985] | 5.196 | 1.98 | 3.596 | 0.5 | no |
6 | [-5.283,6.839] | 10.391 | 3.596 | 4.0 | 0.5 | no |