1 {
3   "name": ""
4  },
5  "nbformat": 3,
6  "nbformat_minor": 0,
7  "worksheets": [
8   {
9    "cells": [
10     {
11      "cell_type": "markdown",
13      "source": [
14       "Described in Goldman, used as argument that Beziers are fractals because they are fixed points of an iterated function system.\n",
15       "\n",
16       "The de Casteljau algorithm splits a Bezier curve into two Bezier segments.\n",
17       "\n"
18      ]
19     },
20     {
22      "level": 3,
24      "source": [
25       "Bezier Curve Definition"
26      ]
27     },
28     {
29      "cell_type": "markdown",
31      "source": [
32       "\n",
33       "Bezier curve $P(t)$ of degree $n$\n",
34       "\n",
35       "$P(t) = \\sum_{j=0}^{n} B_j^n (t) P_j \\quad \\quad \\quad 0 \\leq t \\leq 1$\n",
36       "\n",
37       "For control points $P_0 ... P_n$, with Bernstein basis functions $B_j^n(t) = \\left(^n_j\\right)t^j(1 - t)^{n-j} \\quad \\quad j=0,...n$\n",
38       "\n",
39       "Binomial Coefficients: $\\left(^n_j\\right) = \\frac{n!}{n!(n-j)!}$"
40      ]
41     },
42     {
44      "level": 2,
46      "source": [
47       "Helpers"
48      ]
49     },
50     {
51      "cell_type": "code",
52      "collapsed": false,
53      "input": [
54       "from operator import mul    # or mul=lambda x,y:x*y\n",
55       "from fractions import Fraction\n",
56       "\n",
57       "def nCk(n,k): \n",
58       "    \"\"\" Binomial Coefficient nCk \"\"\"\n",
59       "    # Based on:\n",
60       "    # http://stackoverflow.com/questions/3025162/statistics-combinations-in-python\n",
61       "    # http://stackoverflow.com/questions/279561/what-is-the-python-equivalent-of-static-variables-inside-a-function\n",
62       "\n",
63       "    \n",
64       "    # Dynamic programming is possible in python!\n",
65       "    # Not sure if it's actually worth it though\n",
66       "    try:\n",
67       "        return nCk.dynamic[(n,k)]\n",
68       "    except KeyError:\n",
69       "        nCk.dynamic.update({(n,k) : int( reduce(mul, (Fraction(n-i, i+1) for i in range(k)), 1) )})\n",
70       "    return nCk.dynamic[(n,k)]\n",
71       "nCk.dynamic = {}\n",
72       "\n",
73       "# Simpler version...\n",
74       "def nCr(n,r):\n",
75       "    f = math.factorial\n",
76       "    return f(n) / f(r) / f(n-r)"
77      ],
78      "language": "python",
80      "outputs": [],
81      "prompt_number": 216
82     },
83     {
84      "cell_type": "code",
85      "collapsed": false,
86      "input": [
87       "from time import time\n",
88       "t0 = time()\n",
89       "nCk(100,50)\n",
90       "print str(time() - t0)\n",
91       "t0 = time()\n",
92       "nCr(100,50)\n",
93       "print str(time() - t0)"
94      ],
95      "language": "python",
97      "outputs": [
98       {
99        "output_type": "stream",
100        "stream": "stdout",
101        "text": [
102         "0.000294923782349\n",
103         "0.000556945800781\n"
104        ]
105       }
106      ],
107      "prompt_number": 215
108     },
109     {
111      "level": 2,
113      "source": [
114       "Bezier from Definition"
115      ]
116     },
117     {
118      "cell_type": "code",
119      "collapsed": false,
120      "input": [
121       "def B(n,j, t):\n",
122       "    \"\"\" Bernstein Basis Function of degree n \"\"\"\n",
123       "    assert(0 <= j <= n)\n",
124       "    return nCk(n,j) * t**j * (1. - t)**(n-j)\n",
125       "\n",
126       "def Bezier(P, t):\n",
127       "    \"\"\" Apply Bezier definition to produce a point on the Bezier curve with control points P \"\"\"\n",
128       "    n = len(P)-1\n",
129       "    return [sum([B(n, j, t) * P[j][i] for j in xrange(n+1)]) for i in xrange(len(P[0]))]\n",
130       "        \n",
131       "def PlotBezier(P,nPoints=50,style='o-'):\n",
132       "    \"\"\" Plot a Bezier \"\"\"\n",
133       "    points = map(lambda t : Bezier(P,t), linspace(0,1,nPoints))\n",
134       "    x = [p[0] for p in points]\n",
135       "    y = [p[1] for p in points]\n",
136       "    plot(x,y, style)\n",
137       "    plot([p[0] for p in P], [p[1] for p in P],'ro--')\n"
138      ],
139      "language": "python",
141      "outputs": [],
142      "prompt_number": 266
143     },
144     {
145      "cell_type": "code",
146      "collapsed": false,
147      "input": [
148       "title(\"Bezier Curves of Degree 2\")\n",
149       "PlotBezier([(0,0),(0.4,0.8),(1,1)],nPoints=100,style='-')\n",
150       "PlotBezier([(0,0),(0.4,-0.8),(1,1)],nPoints=100,style='-')\n",
151       "legend()"
152      ],
153      "language": "python",
155      "outputs": [
156       {
158        "output_type": "display_data",
160        "text": [
161         "<matplotlib.figure.Figure at 0x3c4ee50>"
162        ]
163       }
164      ],
165      "prompt_number": 267
166     },
167     {
169      "level": 2,
171      "source": [
172       "Fractal Beziers"
173      ]
174     },
175     {
176      "cell_type": "markdown",
178      "source": [
179       "$L P = Q$ and $M P = R$"
180      ]
181     },
182     {
183      "cell_type": "code",
184      "collapsed": false,
185      "input": [
186       "def L(n):\n",
187       "    try:\n",
188       "        return L.dynamic[n]\n",
189       "    except KeyError:\n",
190       "        result = asarray([asarray([ nCk(j,k) / (2.**j) for k in xrange(n+1)]) for j in xrange(n+1)])\n",
191       "        L.dynamic.update({n : result})\n",
192       "        return result\n",
193       "L.dynamic = {}                        \n",
194       "\n",
195       "def M(n):\n",
196       "    try:\n",
197       "        return M.dynamic[n]\n",
198       "    except KeyError:\n",
199       "        result = asarray([asarray([ nCk(n-j,n-k) / (2.**(n-j)) for k in xrange(n+1)]) for j in xrange(n+1)])\n",
200       "        M.dynamic.update({n : result})\n",
201       "        return result       \n",
202       "M.dynamic = {}\n",
203       "\n"
204      ],
205      "language": "python",
207      "outputs": [],
208      "prompt_number": 296
209     },
210     {
211      "cell_type": "code",
212      "collapsed": false,
213      "input": [
214       "points = [(0,0), (0.5,1), (1,0)]"
215      ],
216      "language": "python",
218      "outputs": [],
219      "prompt_number": 298
220     },
221     {
222      "cell_type": "code",
223      "collapsed": false,
224      "input": [
225       "def Casteljau(P, nPoints):\n",
226       "    while len(P) < nPoints:\n",
227       "        Q = dot(L(len(P)-1), P)\n",
228       "        R = dot(M(len(P)-1), P)\n",
229       "        P = list(Q) + list(R)[1:]\n",
230       "    return P\n",
231       "    "
232      ],
233      "language": "python",
235      "outputs": [],
236      "prompt_number": 329
237     },
238     {
239      "cell_type": "code",
240      "collapsed": false,
241      "input": [
242       "points = [(0,0), (0.5,1), (1,0)]\n",
243       "c = Casteljau(points,50)\n",
244       "PlotBezier(points)\n",
245       "plot([p[0] for p in c], [p[1] for p in c], 'go-')"
246      ],
247      "language": "python",
249      "outputs": [
250       {
252        "output_type": "pyout",
253        "prompt_number": 332,
254        "text": [
255         "[<matplotlib.lines.Line2D at 0x5a3a6d0>]"
256        ]
257       },
258       {
260        "output_type": "display_data",
262        "text": [
263         "<matplotlib.figure.Figure at 0x5942790>"
264        ]
265       }
266      ],
267      "prompt_number": 332
268     },
269     {
270      "cell_type": "markdown",
272      "source": [
273       "So, it sort of works but not giving the same curve as the vanilla method"
274      ]
275     },
276     {
277      "cell_type": "code",
278      "collapsed": false,
279      "input": [],
280      "language": "python",