{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "##
Interactive generation of Bézier and B-spline curves.
Python functional programming implementation of the
de Casteljau and Cox-de Boor algorithms
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The aim of this IPython notebook is twofold: \n", "- first to illustrate the interactive generation of Bézier and B-spline curves using the matplotlib backend `nbagg`, and second\n", "- to give a functional programming implementation of the basic algorithms related to these classes of curves. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bézier and B-spline curves are widely used for interactive heuristic design of free-form curves in [Geometric modeling](http://en.wikipedia.org/wiki/Geometric_modeling)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Their properties are usually illustrated through interactive generation using `C` and `OpenGL` or `Java` applets. The new\n", "`matplotlib nbagg` backend enables the interactive generation of Bézier and B-spline curves in an IPython Notebook.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Why functional programming (FP)? " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lately there has been an increasing interest in pure functional programming and an active debate on whether we can do FP in Python or not.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By [Wikipedia](http://en.wikipedia.org/wiki/Functional_programming):" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "> *Functional programming is a programming paradigm, a style of building the structure and elements of computer programs, that treats computation as the evaluation of mathematical functions and avoids changing-state and mutable data. It is a declarative programming paradigm, which means programming is done with expressions. In functional code, the output value of a function depends only on the arguments that are input to the function, so calling a function f twice with the same value for an argument x will produce the same result f(x) each time. Eliminating side effects, i.e. changes in state that do not depend on the function inputs, can make it much easier to understand and predict the behavior of a program, which is one of the key motivations for the development of functional programming.*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Python is a multi-paradigm programming language: it is both imperative and object-oriented programming language and provides a few [constructs](https://docs.python.org/2/howto/functional.html) for a functional programming style, as well.\n", "\n", "[Here](http://stackoverflow.com/questions/1017621/why-isnt-python-very-good-for-functional-programming) is a discussion on why Python is not very good for FP." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Instead of discussing *pros* and *cons* of Python FP we decided to start a small project and implement it as much as possible in a FP style." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before starting let us recall a few characteristics of FP: \n", "- functions are building blocks of FP. They are first class objects, meaning that they are treated like any other objects. Functions can be passed as arguments to other functions (called *higher order functions*) or be returned by another functions.\n", "- FP defines only *pure functions*, i.e. functions without side effects. Such a function acts only on its input to produce output (like mathematical functions). A pure function never interact with the outside world (it does not perform any I/O operation or modify an instance of a class).\n", "- a FP program consists in evaluating expressions, unlike imperative programming where programs are composed of statements which change global state when executed. \n", "- FP avoids looping and uses recursion instead." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this IPython Notebook we try to implement algorithms related to Bézier and B-spline curves, \n", "using recursion, iterators, higher order functions.\n", "\n", "We also define a class to build interactively a curve. The code for class methods will be imperative \n", "as well as for the last function that prepares data for a closed B-spline curve.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Bézier curves. de Casteljau algorithm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "A Bézier curve of degree $n$ in the $2$-dimensional space is a polynomial curve defined by an ordered set $({\\bf b}_0, {\\bf b}_1, \\ldots, {\\bf b}_n)$ of $n+1$ points, called control points. This set is also called control polygon of the Bézier curve.\n", " \n", "\n", "Its parameterization, $p:t\\in [0,1]\\mapsto p(t)\\in\\mathbb{R}^2$, is defined by:\n", " $$p(t)=\\sum_{k=0}^n{\\bf b}_k B^n_k(t), \\quad t\\in[0,1]$$\n", "where $B^n_k(t)=\\binom nk t^k(1-t)^{n-k}, k=0,1,\\ldots, n$, are Bernstein polynomials.\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To compute a point on a Bézier curve theoretically we have to evaluate the above parameterization $p$ at a parameter $t\\in[0,1]$. In Geometric Modeling a more stable algorithm is used instead, namely the *de Casteljau algorithm*.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "De Casteljau algorithm provides a procedural method to compute recursively a point, $p(t)$, on a Bezier curve. \n", "\n", "Given the control points $({\\bf b}_0, {\\bf b}_1, \\ldots, {\\bf b}_n)$, and a parameter $t\\in [0,1]$, one computes in each step $r=\\overline{1,n}$ of the recursion, the points:\n", " $$ {\\bf\n", "b}_{i}^{r}(t)=(1-t)\\,{\\bf b}_{i}^{r-1}(t)+t\\,{\\bf b}_{i+1}^{r-1}(t),\\:\n", " i=\\overline{0,n-r}, $$ \n", "i.e. the $i^{th}$ point ${\\bf\n", "b}_{i}^{r}(t)$, from the step $r$, is a convex combination of the $i^{th}$ and $(i+1)^{th}$ point from the step $r-1$:\n", "\n", "$$\n", "\\begin{array}{lll} {\\bf b}^{r-1}_i&\\stackrel{1-t}{\\rightarrow}&{\\bf\n", "b}^{r}_i\\\\\n", "{\\bf b}^{r-1}_{i+1}&\\stackrel{\\nearrow}{t}&\\end{array}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The control points calculated in the intermediate steps can be displayed in a triangular array:\n", " $$\n", "\\begin{array}{llllll} \n", "{\\bf b}_0^0 & {\\bf b}_{0}^1 & {\\bf b}_{0}^2& \\cdots & {\\bf b}_0^{n-1}& {\\bf b}^n_0\\\\\n", "{\\bf b}_1^0 & {\\bf b}_{1}^1 & {\\bf b}_{1}^2& \\cdots & {\\bf b}_1^{n-1}& \\\\\n", "{\\bf b}_2^0 & {\\bf b}_{2}^1 & {\\bf b}_{2}^2 & \\cdots & & \\\\\n", "\\vdots & \\vdots &\\vdots & & & \\\\\n", "{\\bf b}_{n-2}^0 & {\\bf b}_{n-2}^1& {\\bf b}_{n-2}^2& & & \\\\\n", "{\\bf b}_{n-1}^0 & {\\bf b}_{n-1}^1 & & & & \\\\\n", "{\\bf b}_n^0 & & & & & \\end{array} $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The points ${\\bf b}_i^0$ denote the given control points ${\\bf b}_i$, $i=\\overline{0,n}$.\n", "\n", "The number of points reduces by 1, from a step to the next one, such that at the final step, $r=n$, we get only one point, ${\\bf b}_0^n(t)$, which is the point $p(t)$\n", "on the Bézier curve.\n", "\n", "\n", "\n", "\n", "The image below illustrates the points computed by de Casteljau algorithm for a Bézier curve defined by 4 control points, and a fixed parameter $t\\in[0,1]$:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAFUCAYAAABoRYRBAAAABmJLR0QA/wD/AP+gvaeTAAAACXBI\nWXMAAAsTAAALEwEAmpwYAAAAHXRFWHRTb2Z0d2FyZQBHUEwgR2hvc3RzY3JpcHQgOS4xMJremEEA\nACAASURBVHic7Z1rbGP5ed6fd+e2u3M9s/bu2p69cWB7ZHsTx5ymdhp7bZeyPHZsd9xS2bgOYPQD\nZRVBChQNJLRJUQRpQ37o7UOgiEiCpE4KZ1hUQYBUGYgwnLZxW3SIGkgLNXGHuRmI0yaivc3W9u56\n/v3wf//U4REvh+S5n+cHDGbEOSKPRPLhe573JsYYEEIIiZ+H0j4BQggpC6fTPgFCSPqISB3AAIAH\noG+M6aV8SoWEgktIyRGRCoCqMWZbvz4AsJruWRUTWgqEFBARaYjIkYgYEWnOOLwW+LonIsHbSARQ\ncAkpIMaYNoD1kIdX4jwXcgwFlxBCEoKCSwghCcGkGSElQES2YCsQPNhqhLYxpq//3QFQ9x1eBdBK\n9gzLASNcQopPA0DXGLNtjNnQ2+6JSBUAXAmYiHh6W88YM0jpXAuNsNOMkGKilQYHsGK76rvdA3A0\n5vYRASbRQ0uBkJJhjBmICBAoB6PQxg8FlxDiOs0qsB5vzxjTSfmUCgkFl5Dy0geG1sPAGNPSr++L\nSJc+bvQwaUZI8amobwtgWLEAHFcieBi1F/qwlQokYhjhElJs+gA2ADRFxA2nqQFYd7aB/t0Bhgm1\nCgD6uTHAKgVCyBCdu9BhAi0eKLiEEADDxBlHM8YIPVxCyIjYikjF1eSSaGGES0jJUXG9A61aAFAx\nxlxP8ZQKCwWXEEISgpYCIYQkBAWXEEISgoJLSMkQkaau37kvIo20z6dM0MMlpERoNcKdwM1X2cab\nDIxwCSkX4yLa+pjbSAwwwiWk4GjZV13/PA7gQuAQRrgJwQiXkAIiIlX1au8DuAc7H6EF4Gn9+yU9\ndINimxyMcAkpCCJSwXEkW4UdSNOFnY1wQlRF5AhWcDn7NiE4LYyQHDNGZHuwQrsaInLtwE4Oo+Am\nBCNcQnKGjlB0IlvDsch2fJt4w9xPHcCuMeZqLCdKTkDBJSQH+ES2pn/3YUW2PY/Ijrlf2goJQkuB\nkAyjUWgNtpzLiezNCEco0lZIEEa4hGQMn8i6+lhnF3RjeizaCglBwSUkA/hqZV1jQgdAN4lLfbUV\n1uMQdDIKBZeQlAg0JFRwLLLthM/jDuzw8e0kH7eMsPGBkATRbQpbYxoSrhpj1pMWW6ULtvcmAiNc\nQmJmQq1sGxMaEpJGKyCOEG0yjoyBVQokdkRQg00CDYxBK+3zSQJfGVcDow0J68uUccWBMWYgIh3Y\n86XgxggFlySBi+K8VM8iZqbUymZOZMfQBbAFgD5ujNDDJXMjgoYIjkRgRNCcdbwxxY6aRKSuiacj\nWNHqw16eXzfGbOdAbAH7wcBtvTHDCJfMjTFoi6AP4CDtc0mLCbWyufVAaSskAwWXkJBo9NfAqMgW\nqS22B1YrxAotBUKmMGauLGBF9qoxpkhiC9gPkKpWVZAYYIRLlkYEW7AJMQ82QdY2BnnwLceiguMi\nWdeQ0EJGyrjiwhjTFxEX5ZaimiRpGOGSZWkA6BqDbWOwobfdE8Ew+aL/rgCoaolY5vA1JNwDcB+2\nlGukIaHIYuvD+bgkBtj4QBZChfMAVmxXfbe7IvqR27OIlnG5SNZfKzvXXNkiodH9fQDXy/o7iBNa\nCiRSjMFABAAyG8mOq5VtIx+1srFDWyFeKLikFIhIAycbEnJbxhUzbkYuBTdiKLgkLobRYlqtvWNq\nZdugyIahC6ApIl5JfOvEYNKMLEtFfVsAw4oFYDQ6Sqy1V0RqIrKrM1539WZXxrVNsZ2N/o76YPIs\nchjhkmXoA9gA0BTBAFZQawDWjTle2WIMeiLxvXkDc2U9FK8hIQ2crZDGuMjCwioFkghu5oIx0QxH\nmTS8GwWvlU0K/f3egy2L4+8zImgpkNwQqJUdO7yb4hANtBXigZYCyTQThnfnZeRh3qGtEDEUXJI5\nxtTKUmTToQPgHqsVooOCS2LH19rriaBmDE5sh50yvJtlXClhjOmJyAD2OWECMgKYNCOpEqiVHeC4\ntZYimwFEZBcAjDEbs44ls6HgksSZMLy7Y4w5EfmSdNHnatcYczXtcykCFFwSKSJwow27futgQq1s\nl7Wy2UebSFjXHAEUXBIZImjguLsLwB//W+DpryBQK2uMYdY7R9BWiA4KLokMERxgZErYAMDbfw34\nk32wISG30FaIDgouWZrjWtkv/X3gPZeO/+eVl4GzTxkDCm3Ooa0QDew0IwshIp6INETkAHZgdR34\nciDp9dKvA9LQ7jDuyco3XQBcob4kjHBJaKbUyrZ1cPUW8Bt3gY+8iJNJswZssqzHagTL7bu3nLft\nAejtre1nNnrU52/LGHM97XPJMxRcMhNfGVcDxyI7Uiurx/Rn1c/6qhV02WQ5fd3bd2/VAGBvbb+r\nX98HcHNvbT+Tvw/9sD0CG1GWgp1mIVAxceMHZ4pKEZhQKzv2zaZ2gRfm96LH9Fy0rN9biEaH23dv\nNQA0YV8nrb21/WmT0TzYS3QX7fcDX2cKY8xARNyCydw/V2nBCHcGbmW2MWZbvz4wxmR6OeKiaPTp\nlioCIWplVTi33O9nwcetwYrNIO8lYxq5HmC24Pq/x4OdfpbZCBegrRAFpU2aacLnSESMiDSnHBpc\nhthTgSgEIlIVkaaI3Id903s43pAQJivdwJK7r4wxXWNMC0BXE2xlS7JtAVjPstgqHQAV/WAmC1Ba\nS8EY0xaRPmw0Mo3CvfFd1I7RhoQW5qyVVduhG5UPq5PAWnrfDY2e+0UuRdLEWWdvbT/zl+m0FZan\ntIJbNibMlZ1bZH33VwWGnmzkOGtBh443kcMk2+27t7Zgrxg86Pnvre33ff9fB9DfW9vv3b57qwLA\ny4HwuhXqZAEouAXGV8bVQITDu/V+a2oDxIqe57YvyZaX0rIGgFUnoLfv3moCuHf77q1VFdgqbIKt\nf/vuLQCo7K3t58Eb7cBu9K1wNvH8UHAVW0M6Go3oC8pdQjmqWNKzjJMptbJRDu9eKkm2CBrZuqi3\nps/XANltGe4FotUWrFfbxLEQ50FgR9B6axflZvZ9kFVKmzQLoNOtzLZvQMc9Eam6S2btrKrCRleZ\ne4OLSF1E7sDWSm7BCu1NY8x1/bkiEVsVulQrCXxJtg5s1NvMepLNlxArQsI1GISQkDDCtfQCXuRI\nNGKM2fZ5lolGdtMYUyvbRoyF6Vqd0cvKpWQg6q27eum8lJZpCVkNwGBvbT9P0SJthQWh4I5Bs7GA\nLxrJSmG+ip5LfgH2xR/7UBG1KqpJ+LaL4H7+QJKtk0FB8J+Pi3q9NE5kUWgrLA4thRCoZ9jUy+k0\nHt9fK+vK2OaplY2CrayKrR9jTF8tlBaAmtb0pnEZX9GGBgDDigXAJ1A5qEiYBm2FBWCEOx0XjSQe\niURVKxvRuWwhh5GMr7SsqlFvH8n8/voANgA0b9+95VrCa7DNDUWpKe7C2grc6DsHFFxLxf/C8UWy\nLWC4vTT2T/MxtbJdpCSyvnOKtLkhDZKc36DDaFz1QdZL1xZG3xN92NdqLjzzLEDB9UUjuhJ6GI0k\ncamuIuAi2chqZSM6t1ibG5Imh6VlWacD+16h4Iak1IKrxfOJRyNjamUzI7KOJJsb0kCf+66zbjRJ\nmsUkW5bpwJZP0lYICZNmCaLzAcbVyt40xrQy9mZvoASRiybZWvrB4pKTiSSDtNusAqDq5uPmCb3y\ncbYCCQHHM4bETRSbtw43UCs7wJjh3VlDx/AFa5NLg89LBxj1TkXfFxVjzHra55IHKLghUC/Ttf62\nZvXxTxjePXWubFbQn7WSh3NNgkmrgXTYTANA121tKCP6erkH4CpthdlQcCPCtzqmDvsGzY3IOlwC\nr6i+7TL4nt/+rc9/8NFzV87+c99/t/fW9jcmfGvh4Ubf8FBwlyAgsq5WtoucZr1FpJml1uUsIiLe\nh37phd9+9IlHVnw3DwBcz8EA8VgQkV0A8M0hIRNg0mxOtHV0S0TuwV5KVWBrZa8aY9aNMbma2erI\na3ND0hhjBqfOnfqm/7YHrz44m9b5ZIQumDgLBQU3gHYl7YrIHZetDojsfdgXVwfA9TyLrEN/zkxO\nQcsSWmWy9fX/9dJv+W9/7Zuv/fqvffg3GyVcDQRgZI4FRXcGtBQCqB/lb+HtwibAXK1sobLWKhC1\nvEzYShpfzXQFxzOSJybNfEm2Qq8GCkJbIRwUXB/6ZtkN3PxlAH+9SCLrh77teHwdgMACq30CpWW5\nvgIKg0a3u8aYq2mfS5ah4PrQN8n9wM0bRY3+3DDxoovBPPiGBvUQwQwJX4Scl9VAC8NqhdlQcAPo\npZGLbHqwA8gLJ0gakfTL2twQREc4VoDjKWMxPUYVBZ3foO+dAa+YJkPBDaARyRGATxtjfiXt84kD\nNjccox88FSQYfU7yhfOOWnJbxpjc7WpLCgpugKJ7UfpmT3wJZNZQcUhd8HyCn5vVQJPwBSuxrXnK\nO6WeFjYBN4e2qJRWbAOJsExUmwRWA7k5zJk4t3nR1VRuEwQFdwyMcAOIyAFssqRwTQAa1XXz+GZe\nBl8iLKmND0uR59Iy2grToeAGEBGDAl4SacLGy9sbeBn8Sao8Xq775zcgBx8UAG2FWVBwfegb9E7R\n/NuyDaXR57EGG83n3h4KJNkyPdoTAHTmc7+s1tU0KLg+dLZn1Rizmva5RElZmhuykgiLkzyUlunz\n0DDG3Ez7XLIGBdeH+re9IolT0ZsbfImwAUrkT2cxAejwNRBdz9J5ZQEKrg/1b1eLcBkKFLu5oWyt\ns9PQ57kK+1xnwqvWQU+dsthYYaHgKm5yvTFG0j6XKNCfp5qVN2BU+BNJRfvZliVLq4H0yqpOW2EU\nCq6iL5BaEfzbIibJipYIi5tJq4ESfHzaCmOg4CpFyqwWybf1JcIyn53PImmWltFWOAkFVynKpKMi\n+LaBCVuZSgjllTRKy2grnISCi+JsHs37UBo3DB1WaOeO0AcrmxVYQfH/mSTWYy+zvcOdwtsVPnsm\nVh/cZyvk+n0VJRRcFKNuMM++rb/1dpYADFY2nSBX9aYqjsUVsOVhPf3bL7YVjG7ycFQn3B6kq/fb\nA9D1DndyLyBJJNlE5D6AFhOcFgouijHHM4/NDeMSYYOVTSemTgidsNZ83+oE1V0Wd4Foo9PBymZQ\niN05OcHv+f94hzu5tXCA+OY3aDNRxRizHtV95hkKLoafwts5vhTPVZIsmEFXK8DVktZxLKZ9HEeq\nfVhhS/1n1PN1A8v9nV9d2PPs5tWaiLq+OS67TgTuw3pgTH62TZdecPPuM2mUiKyXSgWHbh/d+Kz7\n2glWD1awOnmLFn1RuftZXHTsfiYXBecm+RflaqA4bAWR4YczjEFuruwouDkeJ5eHjbv+FtTffm79\nT1fOXX0PjqPDDnQbcp7EKAxqSfijYJfA6wJo5+lDZdnSsrhsBRE0AQpurtAXg2dyuN45y76tiFTe\ncPr8j73//LULP/OGD5zDsVXQgRWdQiSewqJRsLsMbuA4ms/sh2WQRVcDhbEVRNAA0ISNqFthRJSC\nm0O0OLud5ShxHPpB0cqaDbJ++S2bD2A+9sL5a2/99OUbLqrrIIdWQVyoB9zA8dCdDmzUm5sof97V\nQGHyJOrLHoCCW0x8w5Jz1X6YteaGwcqm9+P/+0v/+qkzF6s3znqvf+H8tS6OI7jc/F7TYLCy6YS3\nCqAN+zvLtB/vJ2xpmVYCYdqVZBkEt+w7zVzxd25EQS/PkLTY+qKyKjRi3Xvp/jP//dt/9i9fMQ++\n929evtFbOXf1x2EFI1NRd5ZRS6Gtnm8DwMFgZbMHG/Fm/qpL3zstwOZDNIgZV1rWBbALIHfWXZSU\nPcLNVY1gms0Ng5XNOziOZPCFl//4j37nW3/29Keu3Pji60898o/zFJVlGfV6t3BcIdBG/uyGEzvk\nRFAB/un/AF78PeBNP2Nvx2D0+44jXFirxdM/AwBtY44bWbRKYQvHnm8uXn9lF9xcDddIy7fV6Pa+\n/7YHMN98CPJxCm18qN3gSuecz5ub3/doku1rVeAJf/PKdrB+1ie4AwCrxtjGFrUOGv7b8krZBTc3\nCyPT3LjbuPr8j/3E6//yT1146MxZ381t73Cn1JeHSaEfeC7qHQBo5cFucNjodvQDG0DfGFwPHOcE\nt2sMVn23u1zLyO15pLQerjYMDJYV29t3b7lsrQegt7e2H3m3mu9cExXbv3H5zT/5xtMXXvzIhWff\nfPGhM//FADcAXMZx5QFJALUTNgYrm9uwkV5To9+NPFR+GIO+CPo4nncBTBgeNOH7B2LXAtRmHJp5\nShvhRuHf3r57qwYAe2v7Xf36PoCbe2v7EbYwJu/bfurKjc9deujsDzS8d8hzZy//LNRDdB1Vebqs\nLSL6PLjL7JZ3uJP5LL0ItjD0XP/868Bj9aDvOinC1f8zAGAMRL8eCXSMyUcAUNoIF8ftpCe4fffW\nSBH23tr+pBe0a+l0L5x+4Oso2EqiuUFEvI9efO7n3nHusRd+6vHvO+WdOrcL+2Yefnjovym2KaPP\nw8ZgZbMDYHewslmHjXYz+9wYg5YI2sA/+Angn3xiwVbhPjAU5uEMBRHcF0E3mITLImWOcKcujNTo\n9QDTBdd/vAfbTRNZhJtEkkxEKqsXnv75D55/6u31S29++eqph38VAaEl2cVX1bAFm93P9HPnq30/\nkTvxRbh9+/9WQDU6bgLYMAZtjW6rrv5WZFi7m9kPHEcpI9yYBr5sAViPUGzriHEliojU3nHusY//\n5OPv+fiPXP3uSwC2F0nEaHQ1gPXXet7hTi4u7YqCiuv2YGWzCytK9cHKZmajXVsiJh3YBOC4K8w+\nbK1uU2RYGlYDsO5sA/27AwwTapUJ95U5Sim4iPiyXxNnnb21/UiedK1j9OKonnArtf/Zk+975jNX\n3vYi7At3Y5GoSIv1697hzjqA7mBl8whMpqWCCuzNwcpmE7Z5YuHnNQG6sAHKyJWjRqjXfceEYQtW\njLP4c56glJZCmIWRfksBY4qw99b2+3pcHUB/b22/d/vurQoAbxnhjStJ5pYx/vQTf+VLG97z/xD2\nZzkRCWn2e+hfz0rIDFY2q97hTk9Llw68w53cTV0rGvpB2IQNLDaydtUxzVaY736gLe75iG6B8gru\nEYD1aZaCT3AHAFadiN6+e2tYhK2H3sHxKpfK3tr+UoIT5TBx32jEAYDu0Y3PNhDC69M1NgcImQHX\n42vggJpMMVjZdN5uF9Yyyky32rJbsv1iq3W+Xh6Et3SWgs4i8Obwb3uBiLUF+yJu7q3trwKILKLT\ny/1uBFP23UCRAeyw7yqsgALAatT+nne40x2sbA5gxZ3NEBnBO9xpuUoGAPcGK5st73AnK12VY22F\nMGhbbxNAX+tzK8EmiqzyUNonkAI1LGGw+5JikRZh+z4IlrEjqhoh14wxraMbn+0c3fjsFqzYdgDc\njHjvl6dWAjSyrWi0SzKCd7jT9w53VmGFbWuwsnngnrOU6QKoaHAwF8agZwyuG4NV/ZMLsQVKGOEi\nnoSZ83j7i/i3eulfW9S3DSxjbAHDy/xdPWThqFYvS0f8a9+l6Zb+7aKUMNtvSQp4hzttX7R7X6Pd\n1BomjDF9EenBXollJeqOnbIKbhRJhL4myaquTvf23VsHwEK93gs1N7hEGGz52DZwsi5zyTdVA1as\ne3rfTdhLU3dbC0BNI6YqbFlYJsuRyLCEbF1L+Xb1Q3k1xUoGVx5WGsEtlaWgly8VzBfhVrSpAQBw\n++4tF9W1cNJW6Ll23znOySWxwh7viciWfl/XGLPtbAh9A93T87oZQQQTXP/dgo1im4B9A2sGvAKg\nz2E2+UCfM3cZfqAf0mnQAVBdxFbIK2WLcGuwG0jDfqIPi7Bv3701UoS9t7bf0YqFhVErINT5uIWR\neg4jVQz6htmFRgtxXSp6hzuDwcomEPigYVSbP/S5XIWtsrk3WNlcT7rCpIy2QtkEd+L8hCA6kGbe\nIuzQqIBWZ/m2/mHO43ZHuctDaDtk0m8aX6eZBxvlZr40h1jUSlgdrGzuwka6qyk8f6WyFUplKWDJ\nCoWImdrcICI130Sz7QliuwsbobS9w50kxbavj1+BTg/Ty9SlIn6SDmoFdWBFtz7r+Igpla1QGsHV\nSoB5/dtZBJNvoSLoab6t7oXagp1/uz2uXljLsQ5gI4MovNpJVPz+nlYsAMfnfsLDZllYPlHR3QZw\nR7sNE0FnPPdRgFm3YSiTpRD5wkht561rUq0C2yQx1Y8d19wwuooE7WnnqFHlHf3yeowZ5qF/rU0N\nQ//a1ypaiqikLGjpGGArGNyCyyTowL62crPFYlHKJLhRz6kFAOyt7W/fvnur6v497djgxt3AiumZ\n7bzaI38A+3PENphEk2Cx+dcku6joDmBFt5pQ5UkHwD0R8ZLe15c0ZRLc2D5BpzU7aK1u4y+++vJ/\nBvBWY0wrkAgLlSzwDZVp52HCP8kv3uFOZ7Cy2Yf1dBG36BpjeiLShw0+Ch3llsLD1Uv20BUKUaGb\nI+4D2Lpw7fzeB3/2+z+o/mx1UiJsHNpwsAs7gCRLYruQh02yjyZgV2EbW5Ko1XW2QqEpxbQwrXe9\nY4y5muTjaufZ8EVkvmP+r5ySZ+YZUq6VCHVY7zRzl/f6YdCC9XPrGftAIEuiQusffBTXQPwqbNPO\n1SLbCqWIcKFzBtI+CTklr4U9VisRXNdY5BO+okIFtuL7NykQrlZXvzzQPELkaF7D2QqFpSyCm8ql\n7re//spvBG7qhIlufckxIIVmhnnxDneCLcCkQGgL903Y91BsogsbFMV135mgLIKbSoS7/+IXnoTN\n9rdgh5jPTD74xLanzQyFvbwi+SLQIBGH39pFwSPcwnu46t8eGGMk4cede3ODViK49eS8PCeZRBtg\nmojB6tJtLBvGmEytBYqKMkS4iUe3i2xu8FUibFBsSZbRrREt2K60qJtfCl2tUAbBTdS/nXfjribH\ndmHrcm8m2N1DyMJoUNCBFd0oS8YKbSuURXATiXB9G3fD1te6kpsqfIO+CckJ7kpsd+pRc+CsBL1K\nLByFFtwFFkYuSwMhx8yNqW+k2JJc4TZIwDZHRDkprrC2QqEFFwn6t/P4tkkVkxMSN7rfbhV2QWVU\nU8Z6oODmkkT82+BQmhD4lztSbEmu0auzDejAmwjusgO70bdwNblFF9zYB477Nu6GKmPRBJnzbCm2\npBBosrcFW6O7VOWCXiW6TRCForCC66oFEL+lsDXHxC//XASKLSkUWrnQRTSVC4WsViis4GL+hZFz\n45obwhzrE1smyEiRcd2UyybRCmkrFFlwYy0H823cnblBQpMJFFtSeHyVC/VlKheKaisUWXBj82/d\nfN0w5Wa+dt3E11ATkgZaubCO5SsXCmcrFFJw1b+tIL6EWSjf1ie2G1kdr0hIHOjr3e3EW9QWcLZC\nYXbnFVJwEcPCSMe0jbt+dN20E1u265LSoa97N11s7iSa2go9FCjKLargRr0OHYBdYY4QzQ36iU6x\nJaVHRzq6ObqLVC4UysctquBG7t9qtnQwq7nBN8+2TbElBID1c4HFKhc6AKpFsRUKJ7i+hZGRRbhh\nmxt8YtvhiEVCLFq5sAFbubA1z/eqLVgYW6Fwggsb3Q4i9m8bmFFvq901TmxjXStNSN7QCp112CTa\nvOJZGFuhiIIbdXQ707dVb+oOKLaETEQrF7Yx/8yFwtgKRRXcSPzbML6tf/IXxZaQ6ei2iLkGl/ts\nhdxPECui4EYyknGOoTTDMYvLPiYhZUADkwGAefzcLii42ULbbWdWEoRkZnODzkfwwMlfhMzLNmwn\nWlgR7QCoayCUWwoluIguup3Z3KBdZA1w8hchc6N+bgshS8U0iOoj58mzogluFfZJWRjd3DB1ypga\n/k0A25yPQMjCtAB4cwy5yf3qHTHGpH0OkSEiBsDqojvMNAtam7YE0pck63uHO+uTjiOEzEZLxO4A\nuK5DbyaiSex7AK7GOXY1TgoT4frW3CxjKYTZuNuE9W1ZkUDIkniHOx3YyHXm5t8i2AqFEVws6d/O\n4dtyYwMh0bIBoBpylGMH1jrMJUUS3IXrb8Ns3KVvS0g8aPDSgu1Cm1WFkOs23yIJ7kIRbpiNu/oi\n2IXtJONAGkIiRhsiephhLbi6eA2SckchBFdF08OcEa7W9NVDNDe4LCoH0hASH27AzaxKhNxWKxRC\ncLH4wsgtY8xUEfX5thv0bQmJD61SaMHOWphmLeR29U5RBHfugeM6lGbWBDD6toQkiG+s6cS23zzb\nCkUR3LkGjmsL8NQVPPRtCUmNDdi232nVCLm0FXIvuL6FkaEi3Dk27tK3JSQFfG2/0xJouRxmk3vB\nxfHCyLD+apihNPRtCUmXFoDKlLbfLuxG31zV5BZBcEMPHA/Z3EDflpCU8a3laeg2lRE0wMrdJoii\nCO5MYQzZ3EDflpCMoG2/XUy2FnJXrZBrwQ27MNLV6YaYk0vflpBssY3Jbb8d5MxWyLXg4ti/nVht\n4NvcMKsEjL4tIRnDV5t7ou03j7ZC3gU3jH8bduMufVtCMsiMtt9c2Qp5F9wapgwcD7NxV2kC6NG3\nJSSzbGN826+zFXKx0Te3gjvLvw2zcRcA9Amsg/NtCckseuXZQqADTYOp3EwQy63gwortWEENu3HX\nV5XQmjVtnhCSOi0AtQlRLgU3ZqaNY5zZ3OCOA0b6twkhGUWT2W2cnLPQAVDNg62QZ8EdW38bprkB\nGCbKtkArgZA8cSLK1SqlXNgKeRbcEwNrdChN2DGNrsFh6bXqhJBkUOtvUpRLwQ2LCGoiaIpMHst2\nfKzUgNGFkW6ITZglkroptAo2OBCSR8Z5ubmwFTIjuABcVDprpxEw3r8Ns3HXJcqaANpMlBGSP8ZF\nuXmxFWITXBE0RHAkAiOCSRN/hhgz13qcCnx2gog0EcK3VZgoIyT/jItyu8j4Rt/YBNcYtAGsx3T3\nwwhXh9J0wvi2OgmMiTJCcs4EL7cDoK5loZkkS5ZCKHyDaLphNu4GaIKJMkKKwkiUbjxTdwAADV1J\nREFUqzrQR4ZthdwJLjS6Ddvc4NDhNFUwuiWkEEyJcjO7CeJ0Ug+k1Qee/hkAaBszeQ7CFFz97cyN\nuw5foqzFSWCEFIoWgPuDlc2aXrl2ANwTEW+BLd6xk1SEq0NksG3MMMK8J3JscOu/KwCqIic/oTQJ\ndx/4xieBz69hxgSwAFsA+jp1iBBSEIJRbtZtBTHGxHfnVjgPYMV21Xe7B+AoePuU+3HH+1k3BjPt\nBE2U3QOwSu+WkOKhXaP3oe9xrVqqGGPiStovTCoerjHDmtuwXsu4T6uw3+tqbim2hBQQX5TrtkL0\nkFEfNy9Js3GRbJiOMpcoY80tIcWmBTsvt+IS6VoyminSFtxh0kwEdW3vrfu9XWAYEW8Af/Qq8MrL\nAFqz7AQmyggpD2MqFjJZrZCU4FbUhwUwrFgAtDtMxCbLjEFXhXRMZ5p0gWfOAOe+y5hQEavb4sBE\nGSHloIXjteqZXL2ThOD2YWtfmzqcZhfWa1nXbjTg5CdRb0ylQhUzFkY6tBC6AVoJhJQGf5SbVVsh\n1jpcY9AFcF2/nOa5hpnwE2ZhpMNZCVwISUi5cHW5LRzbCqGao5IgbQ93Hk7Mvx2HJsoqCD/MhhBS\nEAJebuZshVwI7qyFkQG2YMvAmCgjpJx0ADR+9dpHvgLAczNXskBWBDcY8gfX59RgF0ZO9W81uvXA\n6JaQ0qI1993VC09vImObIDIhuG4Wrgg8LQnr+ZojAEa3hJD5aAFofPD8U18GBfckWupV8f3bz9iF\nkX4Y3RJCHC7K/dyb1t4MoJIVWyEzggvYSHfC5odpK9EddTC6JYQc03rkodOfef7h1/07ZCTKzZTg\njkMXRg6mDRnXutsqGN0SQhQX5f6LJ993ARTc0IQpB6N3SwgZR+t7Hn78fVdOncuErZAHwZ3q32p0\nW8N883EJISVAo9z+py/f+ENkYLZCHgR3ln/roluuPCeEjKPzt69+1ylkwFbItOD6lkSOFVxfdEvv\nlhAyifaTp89fe/bMpaqIhBkjEBuZFlwwuiWELInqQ++Hr6x8FSlHuVkX3In+LaNbQsgctP/Wlbc9\nAgruVKZVKDTA6JYQEo7O5VPnHnvnw69P1VbIrOCqf+thjKWgA4brYHRLCAmBlox2Pn3lxteQYpSb\nWcEF3EyFsbvl6d0SQual+0OX33oFQGr1uFkX3EnRbQOMbgkhc+Ad7rTPyqlXPnHxel1HviZOlgV3\nkn+7BaDL6JYQMi+nIJ//8IVnXkJKtkImBVdNbbcIbgijW0LIknR+8PJbLj1x+tFbaTx4JgUXNrrt\nj/FvXXQbdrcZIYQM8Q53uq+Y73z1wxee/Wtp2ApZFVxGt4SQWDgrp37lBy+/+f8hBVshq4I7zr9l\ndEsIiYL2ux95w6PPP/y6TyT9wJkT3HELIwcrmx4Y3RJCIsA73Om/9OCV31s7/8z7k37szAkujv1b\nfxUCo1tCSGRceujsL3zsUuWsiCRqK2RRcMfNT6jj5GZfQghZlPbz5x47/aELz7yY5INmUXBH/NvB\nymYdtsWXgksIiQTvcGfQf+Ub/+m9j77xI0k+bhYFN9hhVgPQ4focQkiUVM5e/sWPX6ycSdJWyJTg\nBhdG+pJljG4JIVHTeezUI6/9yNXvbkR5pyJSF5Ga/j0ytyFTgouTA8frAPpMlhFCosY73Bn8/qvf\n+A9vPee9N6r71C7ZqjGma4zpAGj6/z9rghtMmDXA5ZCEkJh4+7nH2h+98Ny591946v2TjhGRhogc\niYgRkeak45TgosqeXrkDyJ7gDhNmg5XNKqwA004ghMSCd7jTecU8+IsPnX96a9Ixxpg2gPWQdzl1\nuHlmBNd9CvgWRtZhk2WcCkYIiY3ffeXoP1bOXv7+JB4rM4KLk9UJDUxfIEkIIUtz8aGzrbULz1yY\nZitERdYE19kJdcAODE71jAghhecDv/9vvvhbL3/1G6+ZB7/k91snISJbItIUkV39228jBC3QkbzU\n6YjOOQpqOPZJ6mCyjBCSACLSAHBZ/xyISNsYszHh8AaAVVe6qkm0eyKyaozpGWN6Wg7mwfq5I2vC\nMhHhuoWRxpiub0Ekk2WEkCQINj5MW8HTc2KrtGA7YYfVC8aYbWjyTP89JCsRrr+dtw6g5x3uTFqP\nTgghmcAYMxARIFAOFhDlIZmIcDGaMGNnGSEkSYJ605mwLTwUvk6zZrBtOCsRbhVAR2tvK6B/SwhJ\nCGNMW0S60MooX2nqPPSBoT1aN8asA+iKyBF8gp664AYWRjbBQTWEkITR+dvbMw8EKiLiuQhYRFzD\nREvvpyciLf2/CoARLUtdcKH+rTFmoOVgk7KDhBCSJn1YfWqKyAA2WVYDsK5zEwAMRbeG0corAIAY\nYxI835OIyC4AHN34bA/Alne4cz3VEyKEkAhQe6HhLzHLQtLMVSiwFIwQkmtExHONEFqpUMnM8BpX\nHPyPHn/378IKL5NlhJA8swWbfHOM1POm7eHWAPR/9Oo7b8HW3nJQDSEkz7QA1NxcXNj81LDqIW1L\nwdXfspWXEJJ7jDEDTaBVYLePjxQBpB7h/t3H3vXvwSWRhJACMamWN7UIV/3b6meuvO0ZsPaWEFIC\n0oxwa96pc1+/dubCJwGspngehBCSCGkKbrXhPX8fgMclkYSQMpBm0qz60YvPnQG9W0JISUjVUlg5\nd/UlcI0OIaQkpBLhikjthfPXcApyiXYCIaQspBXh1j558fofAPivKT0+IYQkTloebvUD5596GL7l\naoQQUnRSiXC9U+duXjtz4QqYMCOElIjEI1wRqX784vUrAPqcnUAIKRNpWAq1D5y/9lUwuiWElIw0\nBLe6euHp14H+LSGkZCQuuO96+PG1R+T0w2D9LSGkZCQquCJS+djF5y69ah58kcNqCCFlI+kIt/Z9\nj77xpTPy0G8m/LiEEJI6iQruG06ff89feuSJS6CdQAgpIYkK7vvOv+nWyw9e/VPvcIcJM0JI6UhM\ncEWk8u5HnnziW+Y7X0jqMQkhJEskGeHW/ur5p1997NTDewk+JiGEZIbEBPcHLj73nmtnLpwB/VtC\nSElZeJaCiNQBDGAXQPaNMVN92e95+PHa/3ntm3/4lq/8IsvBCCGlZKEI1+1cN8Z0dSVwc8bx3rse\nfv2178DsL/J4hBBSBBa1FGqBr3siErxt5PgXzl/Dk6cf/bkFH48QQnLPUHBFpCEiRyJiRGRqxAqg\nMs+D/L3Hqp96+cGrr7AcjBBSZoaCa4xpA1iP40GeOnPhnb/z7T//b3HcNyGE5IVEqhTee/5Nzz7E\n6gRCSMlZVHCDs2yrmDBu8YeuvPXTz565hO995MlfWPCxCCGkEEwsCxORLdiSLw+2/KttjOkDgDGm\nJyJ1EfFg/dyeMWYQ+H4PQP3ZM5d++kf/5Ivf+uVv/E8Y7MT2gxBCSNYRY8zxF7bS4ABWYFddba0m\n0Rr+2/T2KmAF+MQdW8H2J986xphYPGJCCMkDkwS3a4xZ9d3uATgK3j71jkXu42Q1w3UXJRNCSNkI\n5eH67IJptbZBgkmyPsWWEFJm4mzt7cAKdEWPay/6WIQQUgTmFdw+MNLau61fHwAYsRqMMV0A19Wm\nOJFUI4SQsjHJUqiobwtgmAADgJb+Hbq1V+ctUGwJIaVnXITbB7ABoCkizjKoAVjXQTXAnK29hBBC\nAoLrbAD9kp1hhBASIUlv7SWEkNISe2svIYQQy0KC6+tA87TbjFUIhBAyg5FOs7m/eUprLyGEkFGW\nElxCCCHhYdKMEEISgoJLCCEJQcElhJCEoOASQkhCUHAJISQhKLiEEJIQFFxCCEkICi4hhCQEBZcQ\nQhKCgksIIQlBwSWEkISg4BJCSEJQcAkhJCEouIQQkhAUXEIISQgKLiGEJAQFlxBCEoKCSwghCUHB\nJYSQhKDgEkJIQlBwCSEkISi4hBCSEBRcQghJCAouIYQkBAWXEEISgoJLCCEJQcElhJCEoOASQkhC\nUHAjQkTqIrKV9nnkGf4O84WIbIlIPe3zyBNijEn7HAqBiNwHUAHwwwC+lvLp5JV/BeAN4O8wDzwJ\n4HMA/sAY81zaJ5MXTqd9AgViG0DFGPPLaZ9IXhGRvwP+DnODiLwJQD/t88gTjHAJISQh6OESQkhC\nUHAJISQhKLgRICKeZtjrIuKlfT55RkSaIlJN+zxIOPhczQc93CXRF1zdGLOtYrsFoGWMGaR8arlE\nRAyA68YYJmMyiohUANRhq3IGxpjtlE8pNzDCXZ4tAC0AUJHtAGimekY5RURqAPoU2+xjjGkBYFAx\nJxTcJdCItuqPZo0xPdhPfxISEamJSBP2w8vZCrWUT4tMgB+Ii8M63OWoYnwdoiciVRVfMgNjTBdA\nV0QOYO2YdtrnREgcMMKNhy4AJs/mpwb7uyOkkFBwSSagf0vKAAU3Hhjdzg+jW1J4KLjxMMnbJZOp\nAugBtuxIS48IKRQU3OXowQpFkAEvjefGH+HW+fsjRYSCuwSu7tYfjel80E56Z5Vb2gCqItIArQVS\nUNhptiRai9uEHc/oOs222Wk2Pyylywea4KzBXt15sB+QfZbzzYaCGxH6IhxQMAghk6DgEkJIQvx/\nrWSGQST18YwAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.display import Image\n", "Image(filename='Imag/Decast4p.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We define two functions that implement this recursive scheme. In both cases \n", "we consider 2D points as tuples of float numbers and control polygons as lists of tuples." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First an imperative programming implementation of the de Casteljau algorithm:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "class InvalidInputError(Exception):\n", " pass" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "def deCasteljauImp(b,t): \n", " N=len(b) \n", " if(N<2):\n", " raise InvalidInputError(\"The control polygon must have at least two points\")\n", " a=np.copy(b) #shallow copy of the list of control points and its conversion to a numpy.array \n", " for r in range(1,N): \n", " a[:N-r,:]=(1-t)*a[:N-r,:]+t*a[1:N-r+1,:]# convex combinations in step r \n", " return a[0,:]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "This is a typical imperative programming code: assignment statements, `for` looping.\n", "Each call of this function copies the control points into a numpy array `a`. \n", "The convex combinations `a[i,:] =(1-t) *a[i,:]+t*a[i-1,:]`, $i=\\overline{1,n-r}$, that must be calculated in each step $r$ \n", "are vectorized, computing convex combinations of points from two slices in the `numpy.array` `a`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A discrete version of a Bezier curve, of `nr` points is computed by this function:\n", " " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "def BezierCv(b, nr=200):# compute nr points on the Bezier curve of control points in list bom\n", " t=np.linspace(0, 1, nr)\n", " return [deCasteljauImp(b, t[k]) for k in range(nr)] " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For a FP implementation of de Casteljau algorithm we use standard Python." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we define functions that return an [affine](http://en.wikipedia.org/wiki/Affine_combination)/ [convex](http://en.wikipedia.org/wiki/Convex_combination) combination of two numbers, two 2D points, respectively of each pair of consecutive points in a list of 2D control points:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "cvx=lambda x, y, t: (1-t) * x + t * y # affine/convex combination of two numbers x, y\n", "\n", "cvxP=lambda (P, Q, t): (cvx(P[0], Q[0], t), cvx(P[1], Q[1], t))# affine/cvx comb of two points P,Q\n", "\n", " \n", "def cvxCtrlP(ctrl, t):# affine/cvx combination of each two consecutive points in a list ctrl \n", " # with the same coefficient t \n", " return map(cvxP, zip(ctrl[:-1], ctrl[1:], [t]*(len(ctrl)-1))) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The recursive function implementing de Casteljau scheme:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def deCasteljauF(b, t):\n", " # de Casteljau scheme - computes the point p(t) on the Bezier curve of control polygon b\n", " \n", " if len(b)>1:\n", " return deCasteljauF(cvxCtrlP( b, t), t)\n", " else: \n", " return b[0]\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A Bézier curve of control polygon ${\\bf b}$ is discretized calling `deCasteljauF` function for each\n", "parameter $t_j=j/(nr-1)$, $j=\\overline{0,nr-1}$, where $nr$ is the number of points to be calculated:\n" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "def BezCurve(b, nr=200):\n", " #computes nr points on the Bezier curve of control points b \n", " \n", " return map(lambda s: deCasteljauF(b,s), map(lambda j: j*1.0/(nr-1), xrange(nr)))\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`map(lambda s: deCasteljauF(b,s), map(lambda j: j*1.0/(nr-1), xrange(nr))` means that the function\n", " `deCasteljauF` with fixed list of control points, `b`, and the variable parameter `s` is mapped to the list\n", " of parameters $t_j$ defined above. Instead of defining the list of parameters through comprehension `[j*1.0/(nr-1) for j in xrange(nr)]`\n", " we created it calling the higher order function `map`:\n", "\n", "`map(lambda j: j*1.0/(nr-1), xrange(nr))`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obvioulsly the FP versions of de Casteljau algorithm and `BezCurve` are more compact\n", "than the corresponding imperative versions, but they are not quite readable at the first sight.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To choose interactively the control points of a Bézier curve (and later of a B-spline curve) and to plot the curve, we set now the `matplotlib nbagg` backend, and define the class `BuildB`, below:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "%matplotlib notebook" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "def Curve_plot(ctrl, func):\n", " # plot the control polygon and the corresponding curve discretized by the function func\n", " xc, yc = zip(*func(ctrl))\n", " xd, yd = zip(*ctrl)\n", " plt.plot(xd,yd , 'bo-', xc,yc, 'r') " ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "class BuildB(object): #Build a Bezier/B-spline Curve\n", " \n", " def __init__(self, xlims=[0,1], ylims=[0,1], func=BezCurve):\n", " self.ctrl=[] # list of control points\n", " self.xlims=xlims # limits for x coordinate of control points\n", " self.ylims=ylims # limits for y coordinate\n", " self.func=func # func - function that discretizes a curve defined by the control polygon ctrl\n", " \n", " def callback(self, event): #select control points with left mouse button click\n", " if event.button==1 and event.inaxes:\n", " x,y = event.xdata,event.ydata\n", " self.ctrl.append((x,y)) \n", " plt.plot(x, y, 'bo') \n", " \n", " elif event.button==3: #press right button to plot the curve\n", " Curve_plot(self.ctrl, self.func) \n", " plt.draw() \n", " \n", " else: pass\n", " \n", " def B_ax(self):#define axes lims\n", "\n", " fig = plt.figure(figsize=(6, 5))\n", " ax = fig.add_subplot(111)\n", " ax.set_xlim(self.xlims)\n", " ax.set_ylim(self.ylims)\n", " ax.grid('on') \n", " ax.set_autoscale_on(False)\n", " fig.canvas.mpl_connect('button_press_event', self.callback)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we build the object `C`, set the axes for the curve to be generated, and choose the control points\n", "with the left mouse button click. A right button click generates the corresponding curve:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support.' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " fig.waiting = false;\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " this.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width);\n", " canvas.attr('height', height);\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
')\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('');\n", " button.click(method_name, toolbar_event);\n", " button.mouseover(tooltip, toolbar_mouse_event);\n", " nav_element.append(button);\n", " }\n", "\n", " // Add the status bar.\n", " var status_bar = $('');\n", " nav_element.append(status_bar);\n", " this.message = status_bar[0];\n", "\n", " // Add the close button to the window.\n", " var buttongrp = $('
');\n", " var button = $('');\n", " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", " button.mouseover('Close figure', toolbar_mouse_event);\n", " buttongrp.append(button);\n", " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", " titlebar.prepend(buttongrp);\n", "}\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(el){\n", " // this is important to make the div 'focusable\n", " el.attr('tabindex', 0)\n", " // reach out to IPython and tell the keyboard manager to turn it's self\n", " // off when our div gets focus\n", "\n", " // location in version 3\n", " if (IPython.notebook.keyboard_manager) {\n", " IPython.notebook.keyboard_manager.register_events(el);\n", " }\n", " else {\n", " // location in version 2\n", " IPython.keyboard_manager.register_events(el);\n", " }\n", "\n", "}\n", "\n", "mpl.figure.prototype._key_event_extra = function(event, name) {\n", " var manager = IPython.notebook.keyboard_manager;\n", " if (!manager)\n", " manager = IPython.keyboard_manager;\n", "\n", " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", " event.shiftKey = false;\n", " // Send a \"J\" for go to next cell\n", " event.which = 74;\n", " event.keyCode = 74;\n", " manager.command_mode();\n", " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", "mpl.figure.prototype.handle_save = function(fig, msg) {\n", " fig.ondownload(fig, null);\n", "}\n", "\n", "\n", "mpl.find_output_cell = function(html_output) {\n", " // Return the cell and output element which can be found *uniquely* in the notebook.\n", " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", " // IPython event is triggered only after the cells have been serialised, which for\n", " // our purposes (turning an active figure into a static one), is too late.\n", " var cells = IPython.notebook.get_cells();\n", " var ncells = cells.length;\n", " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", " data = data.data;\n", " }\n", " if (data['text/html'] == html_output) {\n", " return [cell, data, j];\n", " }\n", " }\n", " }\n", " }\n", "}\n", "\n", "// Register the function which deals with the matplotlib target/channel.\n", "// The kernel may be null if the page has been refreshed.\n", "if (IPython.notebook.kernel != null) {\n", " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", "}\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cv=BuildB()\n", "cv.B_ax()" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "left, right=subdivision(cv.ctrl, 0.47)\n", "plot_polygon(left, 'g')\n", "plot_polygon(right, 'm')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Multi-affine de Casteljau algorithm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The above de Casteljau algorithm is the classical one. There is a newer approach to define and study Bézier\n", "curves through polarization.\n", "\n", "Every polynomial curve of degree n, parameterized by $p:[a,b]\\to\\mathbb{R}^2$ defines a [symmetric multiaffine\n", " map](http://www.cs.utah.edu/~xchen/blossom/node7.html) $g:\\underbrace{[a,b]\\times[a,b]\\times[a,b]}_{n}\\to \\mathbb{R}^d$, such that\n", " $p(t)=g(\\underbrace{t,t,\\ldots,t }_{n})$, for every $t\\in[a,b]$.\n", " \n", "$g$ is called the *polar form* of the polynomial curve. An argument $(u_1, u_2, \\ldots, u_n)$ of the polar form $g$ is called *polar bag*. \n", "\n", "If $p(t)$ is the Bernstein parameterization of a Bezier curve of control points \n", "${\\bf b}_0,{\\bf b}_{1}, \\ldots, {\\bf b}_n$, and $g:[0,1]^n\\to\\mathbb{R}^2$ its polar form, then the control points are related to $g$ as follows [Gallier]:\n", "\n", "$${\\bf b}_j=g(\\underbrace{0, 0, \\ldots, 0}_{n-j}, \\underbrace{1,1, \\ldots, 1}_{j}), \\quad j=0,1, \\ldots, n$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This relationship allows to define Bézier curves by a parameterization $P$ defined not only on the interval $[0,1]$ but also on an arbitrary interval $[r,s]$. Namely, given a symmetric multiaffine map $g:[r,s]^n\\to\\mathbb{R}^2$\n", "the associated polynomial curve $P(t)=g(\\underbrace{t,t,\\ldots, t}_{n})$, $t\\in[r,s]$, expressed as a Bézier curve has the control points defined by [Gallier]:\n", "$${\\bf b}_j=g(\\underbrace{r, r, \\ldots, r}_{n-j}, \\underbrace{s,s, \\ldots, s}_{j}), \\quad j=0,1, \\ldots, n$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Given the control points ${\\bf b}_0,{\\bf b}_{1}, \\ldots, {\\bf b}_n$ of a Bézier curve, and a polar bag $(u_1, u_2, \\ldots, u_n)\\in[r,s]^n$ the multi-affine de Casteljau algorithm\n", "evaluates the corresponding multi-affine map, $g$, at this polar bag through a recursive formula similar to the classical\n", "de Casteljau formula:\n", " \n", "$${\\bf b}^k_i= \\displaystyle\\frac{s-u_k}{s-r}\\:{\\bf b}_i^{k-1}+ \\displaystyle\\frac{u_k-r}{s-r}\\:{\\bf b}_{i+1}^{k-1}, \\quad k=\\overline{1,n}, i=\\overline{0,n-k}$$ \n", " \n", "\n", "The point $b_0^n$ computed in the last step is the polar value, $g(u_1,u_2, \\ldots, u_n)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Unlike the classical de Casteljau formula, where in each step the parameter for the convex combination\n", "is the same, here in each step $k$ the parameter involved in convex combinations changes, namely it is \n", "$\\displaystyle\\frac{u_k-r}{s-r}$. \n", "\n", "In order to define a recursive function to implement this scheme\n", "we consider the polar bag as an [iterator](http://nbviewer.ipython.org/github/empet/pytwist/blob/master/Generators-and-Dynamical-Systs.ipynb) associated to the list containing its coordinates, i.e. if `L=[u[0], u[1], ..., u[n-1]]` is the list, `iter(L)` is its iterator:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "def deCasteljauAff(b, r, s, u): #multi-affine de Casteljau algorithm\n", " \n", " # b is the list of control points b_0, b_1, ..., b_n\n", " # [r,s] is a subinterval \n", " # u is the iterator associated to the polar bag [u_0, u_1, u_{n-1}], n=len(b)-1 \n", " \n", " if len(b)>1:\n", " return deCasteljauAff(cvxCtrlP( b, (u.next()-r)/(s-r)), r, s, u)\n", " else: return b[0]\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Usually we should test the concordance between the length (n+1) of the control polygon and the number of elements of the iterator \n", "u. Since a `listiterator` has no length we have to count its elements. Functionally this number would get as:\n", "`len(map(lambda item: item, u))`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What characteristic elements of a Bézier curve can the multi-affine de Casteljau algorithm compute?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**1.** The direction of the tangent to a Bézier curve of degree $n$, at a point $p(t)$, is defined by the vector\n", "$\\overrightarrow{{\\bf b}^{n-1}_0{\\bf b}^{n-1}_1}$. The end points of this vector are the points computed in\n", "the $(n-1)^{th}$\n", "step of the classical de Casteljauscheme. On the other hand, these points are polar values of the corresponding\n", "multiaffine map $g$ [Gallier], namely:\n", "\n", "$${\\bf b}^{n-1}_0=g(\\underbrace{t,t,\\ldots, t}_{n-1}, 0), \\quad {\\bf b}^{n-1}_1=g(\\underbrace{t,t,\\ldots, t}_{n-1}, 1)$$\n", "Thus they can be computed by the function `deCasteljauAff`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us generate a Bézier curve and draw the tangent at the point corresponding to $t=0.45$:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support.' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " fig.waiting = false;\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " this.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width);\n", " canvas.attr('height', height);\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
')\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('');\n", " button.click(method_name, toolbar_event);\n", " button.mouseover(tooltip, toolbar_mouse_event);\n", " nav_element.append(button);\n", " }\n", "\n", " // Add the status bar.\n", " var status_bar = $('');\n", " nav_element.append(status_bar);\n", " this.message = status_bar[0];\n", "\n", " // Add the close button to the window.\n", " var buttongrp = $('
');\n", " var button = $('');\n", " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", " button.mouseover('Close figure', toolbar_mouse_event);\n", " buttongrp.append(button);\n", " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", " titlebar.prepend(buttongrp);\n", "}\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(el){\n", " // this is important to make the div 'focusable\n", " el.attr('tabindex', 0)\n", " // reach out to IPython and tell the keyboard manager to turn it's self\n", " // off when our div gets focus\n", "\n", " // location in version 3\n", " if (IPython.notebook.keyboard_manager) {\n", " IPython.notebook.keyboard_manager.register_events(el);\n", " }\n", " else {\n", " // location in version 2\n", " IPython.keyboard_manager.register_events(el);\n", " }\n", "\n", "}\n", "\n", "mpl.figure.prototype._key_event_extra = function(event, name) {\n", " var manager = IPython.notebook.keyboard_manager;\n", " if (!manager)\n", " manager = IPython.keyboard_manager;\n", "\n", " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", " event.shiftKey = false;\n", " // Send a \"J\" for go to next cell\n", " event.which = 74;\n", " event.keyCode = 74;\n", " manager.command_mode();\n", " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", "mpl.figure.prototype.handle_save = function(fig, msg) {\n", " fig.ondownload(fig, null);\n", "}\n", "\n", "\n", "mpl.find_output_cell = function(html_output) {\n", " // Return the cell and output element which can be found *uniquely* in the notebook.\n", " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", " // IPython event is triggered only after the cells have been serialised, which for\n", " // our purposes (turning an active figure into a static one), is too late.\n", " var cells = IPython.notebook.get_cells();\n", " var ncells = cells.length;\n", " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", " data = data.data;\n", " }\n", " if (data['text/html'] == html_output) {\n", " return [cell, data, j];\n", " }\n", " }\n", " }\n", " }\n", "}\n", "\n", "// Register the function which deals with the matplotlib target/channel.\n", "// The kernel may be null if the page has been refreshed.\n", "if (IPython.notebook.kernel != null) {\n", " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", "}\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "Bez=BuildB()\n", "Bez.B_ax()" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "br=redefineBezier(Bez.ctrl, 0.3, 0.67)\n", "plot_polygon(br, 'g')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**3.** The function `redefineBezier` can also be invoked to compute the left and right subpolygons resulted from a subdivision of a Bézier curve, \n", "at a point corresponding to the paramater $s$. Namely the left subpolygon is returned by\n", "`redefineBezier(b, 0, s)`, whereas the right subpolygon by `redefineBezier(b, s, 1)` " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### B-spline curves" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We give a procedural definition of a B-spline curve of degree $k$, not an analytical one." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following data: \n", "\n", "- an interval $[a,b]$;\n", "\n", "- an integer $k\\geq 1$;\n", "\n", "- a sequence of knots:\n", "\n", " $$u_0 = u_1 = \\cdots = u_k< u_{k+1} \\leq \\cdots \\leq u_{m−k−1} < u_{m−k} = \\cdots = u_{m−1} = u_m$$\n", "\n", " with $u_0 = a, u_m = b$, $m−k > k$, and each knot $u_{k+1}, \\ldots, u_{m-k-1}$ of multiplicity at most $k$;\n", "- the control points ${\\bf d}_0, {\\bf d}_1, \\ldots, {\\bf d}_{m-k-1}$, called de Boor points;\n", "\n", "define a curve $s:[a,b]\\to\\mathbb{R}^2$, such that for each $t$ in an interval $[u_J, u_{J+1}]$, with\n", "$u_J1$, the coefficients $\\omega^r_i$ are computed with the same formula as for\n", "$r=1$, i.e. $\\omega_i^r=\\displaystyle\\frac{t-u_i}{u_i+k-u_i}$ but for the knots in the list\n", "`u[r-1: -r+1]` and $k=k-r+1$ we \n", "define the function `omega` below, instead of `Omegas`:\n" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [], "source": [ "def omega(u, k, t):\n", " # defines the list of coefficients for the convex combinations performed in a step of de Boor algo\n", " \n", " #if (len(u)!=2*k+1 or tu[k+1]):\n", " #raise InvalidInputError('the list u has not the length 2k+1 or t isn't within right interval')\n", " \n", " return map(lambda j: (t-u[j]) / (u[j+k]-u[j]), xrange(1,k+1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We also need a new function that calculates the convex combinations of each pair of points in a list,\n", "with distinct coefficients given in a list `alpha`:" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [], "source": [ "def cvxList(d, alpha):\n", " #len(d)=len(alpha)+1\n", " return map(cvxP, zip(d[:-1], d[1:], alpha))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The recursive Cox-de Boor formula is now implemented in the following way:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "def DeBoor(d, u, k, t):\n", " #len(d) must be (k+1) and len(u)=2*k+1:\n", " # this algorithm evaluates a point c(t) on an arc of B-spline curve\n", " # of degree k, defined by:\n", " # u_0<=u_1<=... u_k');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " fig.waiting = false;\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " this.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width);\n", " canvas.attr('height', height);\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
')\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('