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       "I think the algorithm developed by Goldman might be worse than just applying the Bezier definition...\n",
17       "\n"
18      ]
19     },
20     {
22      "level": 2,
24      "source": [
25       "Binomial Coefficients"
26      ]
27     },
28     {
29      "cell_type": "code",
30      "collapsed": false,
31      "input": [
32       "from operator import mul    # or mul=lambda x,y:x*y\n",
33       "from fractions import Fraction\n",
34       "\n",
35       "def nCk(n,k): \n",
36       "    \"\"\" Binomial Coefficient nCk \"\"\"\n",
37       "    # Based on:\n",
38       "    # http://stackoverflow.com/questions/3025162/statistics-combinations-in-python\n",
39       "    # http://stackoverflow.com/questions/279561/what-is-the-python-equivalent-of-static-variables-inside-a-function\n",
40       "\n",
41       "    \n",
42       "    # Dynamic programming is possible in python!\n",
43       "    # Not sure if it's actually worth it though\n",
44       "    try:\n",
45       "        return nCk.dynamic[(n,k)]\n",
46       "    except KeyError:\n",
47       "        nCk.dynamic.update({(n,k) : int( reduce(mul, (Fraction(n-i, i+1) for i in range(k)), 1) )})\n",
48       "    return nCk.dynamic[(n,k)]\n",
49       "nCk.dynamic = {}\n",
50       "\n",
51       "# Simpler version...\n",
52       "def nCr(n,r):\n",
53       "    f = math.factorial\n",
54       "    return f(n) / f(r) / f(n-r)"
55      ],
56      "language": "python",
58      "outputs": [],
59      "prompt_number": 2
60     },
61     {
62      "cell_type": "code",
63      "collapsed": false,
64      "input": [
65       "from time import time\n",
66       "t0 = time()\n",
67       "nCk(100,50)\n",
68       "print str(time() - t0)\n",
69       "t0 = time()\n",
70       "nCr(100,50)\n",
71       "print str(time() - t0)"
72      ],
73      "language": "python",
75      "outputs": [
76       {
77        "output_type": "stream",
78        "stream": "stdout",
79        "text": [
80         "0.00108003616333\n",
81         "0.000107049942017\n"
82        ]
83       }
84      ],
85      "prompt_number": 3
86     },
87     {
89      "level": 2,
91      "source": [
92       "Bezier from Definition"
93      ]
94     },
95     {
96      "cell_type": "markdown",
98      "source": [
99       "\n",
100       "Bezier curve $P(t)$ of degree $n$\n",
101       "\n",
102       "$P(t) = \\sum_{j=0}^{n} B_j^n (t) P_j \\quad \\quad \\quad 0 \\leq t \\leq 1$\n",
103       "\n",
104       "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",
105       "\n",
106       "Binomial Coefficients: $\\left(^n_j\\right) = \\frac{n!}{n!(n-j)!} = \\text{nCj}$"
107      ]
108     },
109     {
110      "cell_type": "code",
111      "collapsed": false,
112      "input": [
113       "def B(n,j, t):\n",
114       "    \"\"\" Bernstein Basis Function of degree n \"\"\"\n",
115       "    assert(0 <= j <= n)\n",
116       "    return nCk(n,j) * t**j * (1. - t)**(n-j)\n",
117       "\n",
118       "def Bezier(P, t):\n",
119       "    \"\"\" Apply Bezier definition to produce a point on the Bezier curve with control points P \"\"\"\n",
120       "    n = len(P)-1\n",
121       "    return [sum([B(n, j, t) * P[j][i] for j in xrange(n+1)]) for i in xrange(len(P[0]))]\n",
122       "\n",
123       "def BezierCurve(P, nPoints):\n",
124       "    return map(lambda t : Bezier(P,t), linspace(0,1,nPoints))\n",
125       "        \n",
126       "def PlotBezier(P,nPoints=50,style='-'):\n",
127       "    \"\"\" Plot a Bezier \"\"\"\n",
128       "    points = map(lambda t : Bezier(P,t), linspace(0,1,nPoints))\n",
129       "    x = [p[0] for p in points]\n",
130       "    y = [p[1] for p in points]\n",
131       "    plot(x,y, style)\n",
132       "    plot([p[0] for p in P], [p[1] for p in P],'ro--')\n"
133      ],
134      "language": "python",
136      "outputs": [],
137      "prompt_number": 4
138     },
139     {
140      "cell_type": "code",
141      "collapsed": false,
142      "input": [
143       "title(\"Bezier Curves of Degree 2\")\n",
144       "PlotBezier([(0,0),(0.4,0.8),(1,1)],nPoints=100,style='-')\n",
145       "PlotBezier([(0,0),(0.4,-0.8),(1,1)],nPoints=100,style='-')"
146      ],
147      "language": "python",
149      "outputs": [
150       {
151        "output_type": "stream",
152        "stream": "stderr",
153        "text": [
154         "/usr/lib/pymodules/python2.7/matplotlib/axes.py:4486: UserWarning: No labeled objects found. Use label='...' kwarg on individual plots.\n",
155         "  warnings.warn(\"No labeled objects found. \"\n"
156        ]
157       },
158       {
160        "output_type": "display_data",
162        "text": [
163         "<matplotlib.figure.Figure at 0x2bb2bd0>"
164        ]
165       }
166      ],
167      "prompt_number": 5
168     },
169     {
171      "level": 2,
173      "source": [
174       "De Casteljau Algorithm"
175      ]
176     },
177     {
178      "cell_type": "markdown",
180      "source": [
181       "Subdivide the original $n$ control points $\\{P_0 ... P_n\\}$ into $\\{Q_0 ... Q_n\\}$ and $\\{R_0 ... R_n\\}$ (ie: Makes $2n$ points)\n",
182       "\n",
183       "Defining $(n+1)\\times(n+1)$ matrices $L$ and $M$ such that $L P = Q$ and $M P = R$\n",
184       "\n",
185       "Then it can be shown (subdividing at $t = 1/2$)\n",
186       "\n",
187       "$L_{i,j} = \\left(\\frac{\\left(^j_k\\right)}{2^j}\\right)$ and $M_{i,j} = \\left(\\frac{\\left(^{n-j}_{n-k}\\right)}{2^{n-j}}\\right)$"
188      ]
189     },
190     {
191      "cell_type": "markdown",
193      "source": [
194       "**Confusion:** Do $Q$ and $R$ converge to the same curve as applying the definition to produce points along $P(t)$?\n",
195       "\n",
196       "It is implied in Goldman, but it doesn't look like it, unless I made a mistake. See below."
197      ]
198     },
199     {
200      "cell_type": "code",
201      "collapsed": false,
202      "input": [
203       "def L(n):\n",
204       "    \"\"\" Gives L of size (n+1)*(n+1) \"\"\"\n",
205       "    try:\n",
206       "        return L.dynamic[n]\n",
207       "    except KeyError:\n",
208       "        result = asarray([asarray([ nCk(j,k) / (2.**j) for k in xrange(n+1)]) for j in xrange(n+1)])\n",
209       "        L.dynamic.update({n : result})\n",
210       "        return result\n",
211       "L.dynamic = {}                        \n",
212       "\n",
213       "def M(n):\n",
214       "    \"\"\" Gives M of size (n+1)*(n+1) \"\"\"\n",
215       "    try:\n",
216       "        return M.dynamic[n]\n",
217       "    except KeyError:\n",
218       "        result = asarray([asarray([ nCk(n-j,n-k) / (2.**(n-j)) for k in xrange(n+1)]) for j in xrange(n+1)])\n",
219       "        M.dynamic.update({n : result})\n",
220       "        return result       \n",
221       "M.dynamic = {}\n",
222       "\n"
223      ],
224      "language": "python",
226      "outputs": [],
227      "prompt_number": 6
228     },
229     {
230      "cell_type": "code",
231      "collapsed": false,
232      "input": [
233       "L(2)"
234      ],
235      "language": "python",
237      "outputs": [
238       {
240        "output_type": "pyout",
241        "prompt_number": 7,
242        "text": [
243         "array([[ 1.  ,  0.  ,  0.  ],\n",
244         "       [ 0.5 ,  0.5 ,  0.  ],\n",
245         "       [ 0.25,  0.5 ,  0.25]])"
246        ]
247       }
248      ],
249      "prompt_number": 7
250     },
251     {
252      "cell_type": "code",
253      "collapsed": false,
254      "input": [
255       "def Casteljau(P, nPoints):\n",
256       "    \"\"\" Approximates the Bezier curve with control points P using the de Casteljau subdivision algorithm \"\"\"\n",
257       "    while len(P) < nPoints:\n",
258       "        Q = dot(L(len(P)-1), P)\n",
259       "        R = dot(M(len(P)-1), P)\n",
260       "        P = list(Q) + list(R)\n",
261       "    return P\n",
262       "    "
263      ],
264      "language": "python",
266      "outputs": [],
267      "prompt_number": 8
268     },
269     {
270      "cell_type": "code",
271      "collapsed": false,
272      "input": [
273       "points = [(0,0), (0.5,1),(1,0),(1.4,-2)]\n",
274       "c = Casteljau(points,50)\n",
275       "PlotBezier(points)\n",
276       "plot([p[0] for p in c], [p[1] for p in c], 'g-')"
277      ],
278      "language": "python",
280      "outputs": [
281       {
283        "output_type": "pyout",
284        "prompt_number": 9,
285        "text": [
286         "[<matplotlib.lines.Line2D at 0x2af9b10>]"
287        ]
288       },
289       {
291        "output_type": "display_data",
293        "text": [
294         "<matplotlib.figure.Figure at 0x2af9950>"
295        ]
296       }
297      ],
298      "prompt_number": 9
299     },
300     {
301      "cell_type": "markdown",
303      "source": [
304       "So, it sort of works but not giving the same curve? Or maybe I just need *lots* more points and they will converge."
305      ]
306     },
307     {
309      "level": 2,
311      "source": [
312       "Iterated Function Systems"
313      ]
314     },
315     {
316      "cell_type": "markdown",
318      "source": [
319       "Defining $L_p = P^{-1} L P$ and $M_p = P^{-1} M P$, then $P L_p = L P = Q$ and $P M_p = M P = R$\n",
320       "\n",
321       "\n",
322       "$\\{L_p, M_p\\}$ is an iterated function system, $\\implies$ Beziers are fractals.\n",
323       "\n",
324       "Problem: $P$ needs to be invertable; introduce homogeneous coordinates $\\langle$insert hand waving here$\\rangle$ to \"lift\" it and then it all magically works"
325      ]
326     },
327     {
328      "cell_type": "code",
329      "collapsed": false,
330      "input": [
331       "def lift(P):\n",
332       "    \"\"\" Lift P to make it square by introducing homogeneous coordinates\"\"\"\n",
333       "    n = len(P)\n",
334       "    result = []\n",
335       "    for i in xrange(len(P)):\n",
336       "        ones = min(n-len(P[i]), max(1, i-1))\n",
337       "        zeroes = max(n-len(P[i])-ones, 0)\n",
338       "        result += [asarray(list(P[i]) + [1. for _ in xrange(ones)] + [0. for _ in xrange(zeroes)])]\n",
339       "    return asarray(result)\n",
340       "\n",
341       "def lower(P):\n",
342       "    \"\"\" Lower square P back to list of 2D points \"\"\"\n",
343       "    return [row[0:2] for row in P]\n",
344       "\n",
345       "def Lp(P):\n",
346       "    return dot(inv(P), dot(L(len(P)-1), P))\n",
347       "\n",
348       "def Mp(P):\n",
349       "    return dot(inv(P), dot(M(len(P)-1), P)) # Typo in Goldman? Reads \"P^-1 * R * P\" instead of \"P^-1 * M * P\""
350      ],
351      "language": "python",
353      "outputs": [],
354      "prompt_number": 10
355     },
356     {
357      "cell_type": "markdown",
359      "source": [
360       "So, Goldman's method is:\n",
361       "\n",
362       "<ol>\n",
363       " <li> lift the control points\n",
364       "P from the ambient affine space of dimension $d$ to the\n",
365       "ambient affine space of dimension $n$;\n",
366       " </li><li> generate the corresponding higher dimensional Bezier curve using the\n",
367       "iterated function system $\\{L_p , M_p \\}$ ; \n",
368       "</li><li> project the\n",
369       "resulting $n$-dimensional Bezier curve orthogonally back\n",
370       "down to the original dimension $d$.\n",
371       "</ol>\n"
372      ]
373     },
374     {
375      "cell_type": "code",
376      "collapsed": false,
377      "input": [
378       "def Goldman(P, nPoints):\n",
379       "    while len(P) < nPoints:\n",
380       "        Pl = lift(P)\n",
381       "        LP = Lp(Pl)\n",
382       "        MP = Mp(Pl)\n",
383       "        Q = dot(Pl, LP)\n",
384       "        R = dot(Pl, MP)\n",
385       "        P = lower(Q) + lower(R)\n",
386       "    return P\n",
387       "        "
388      ],
389      "language": "python",
391      "outputs": [],
392      "prompt_number": 11
393     },
394     {
395      "cell_type": "code",
396      "collapsed": false,
397      "input": [
398       "points = [(0,0), (0.5,1),(1,0),(1.4,-2)]\n",
399       "c = Casteljau(points,50)\n",
400       "c2 = Goldman(points,50)\n",
401       "PlotBezier(points)\n",
402       "plot([p[0] for p in c], [p[1] for p in c], 'g-')\n",
403       "plot([p[0] for p in c2], [p[1] for p in c2], 'r-')"
404      ],
405      "language": "python",
407      "outputs": [
408       {
410        "output_type": "pyout",
411        "prompt_number": 12,
412        "text": [
413         "[<matplotlib.lines.Line2D at 0x30a8750>]"
414        ]
415       },
416       {
418        "output_type": "display_data",
420        "text": [
421         "<matplotlib.figure.Figure at 0x2bfe7d0>"
422        ]
423       }
424      ],
425      "prompt_number": 12
426     },
427     {
428      "cell_type": "markdown",
430      "source": [
431       "The algorithm works, I guess, but is it useful? It will be a lot slower than de Casteljau, because it does a lot more matrix manipulation stuff, but it gives no more accuracy on each iteration."
432      ]
433     },
434     {
436      "level": 2,
438      "source": [
439       "Time Performance of Casteljau, Goldman, and Basic Bezier creation"
440      ]
441     },
442     {
443      "cell_type": "code",
444      "collapsed": false,
445      "input": [
446       "from time import time\n",
447       "points = [(0,0), (0.5,1),(1,0),(1.4,-2)]\n",
448       "costs = {BezierCurve : [], Casteljau : [], Goldman : []}\n",
449       "stddev = {}\n",
450       "nAverages = 10\n",
451       "depth=8\n",
452       "for method in costs:\n",
453       "\n",
454       "    stddev.update({method : []})\n",
455       "    for n in xrange(1,depth):\n",
456       "        #display(\"depth %d for %s\" % (n, str(method)))\n",
457       "        c = []\n",
458       "        for _ in xrange(nAverages):\n",
459       "            t0 = time()\n",
460       "            p = method(points, nPoints=len(points) * 2.**n)\n",
461       "            c += [time()-t0]\n",
462       "        costs[method] += [mean(c)]\n",
463       "        stddev[method] += [var(c)**0.5]"
464      ],
465      "language": "python",
467      "outputs": [],
468      "prompt_number": 17
469     },
470     {
471      "cell_type": "code",
472      "collapsed": false,
473      "input": [
474       "title(\"Time to produce $n$-dimensional Bezier Curve\")\n",
475       "xlabel(\"$n$\")\n",
476       "ylabel(\"$t$ (s)\")\n",
477       "x = [len(points) * 2.**n for n in xrange(1,depth)]\n",
478       "for method in costs:\n",
479       "    errorbar(x,costs[method],yerr=stddev[method],label =str(method).split(\" \")[1])\n",
480       "legend(loc=\"upper left\")"
481      ],
482      "language": "python",
484      "outputs": [
485       {
487        "output_type": "pyout",
488        "prompt_number": 18,
489        "text": [
490         "<matplotlib.legend.Legend at 0x2f7f790>"
491        ]
492       },
493       {
495        "output_type": "display_data",
497        "text": [
498         "<matplotlib.figure.Figure at 0x3a1be50>"
499        ]
500       }
501      ],
502      "prompt_number": 18
503     },
504     {
506      "level": 2,
508      "source": [
509       "Approximating a Circle using Cubic Beziers"
510      ]
511     },
512     {
513      "cell_type": "code",
514      "collapsed": false,
515      "input": [
516       "k = 4. * (2.**0.5 - 1.) / 3.\n",
517       "points = [(0,1), (k, 1), (1, k), (1,0)]\n",
518       "c = BezierCurve(points, 50) # This appears to be more accurate than using the de Casteljau algorithm :S\n",
519       "c += list(array([-1,1]) * c)\n",
520       "c += list(array([1,-1]) * c)\n",
521       "scatter([p[0] for p in c], [p[1] for p in c])"
522      ],
523      "language": "python",
525      "outputs": [
526       {
528        "output_type": "pyout",
529        "prompt_number": 125,
530        "text": [
531         "<matplotlib.collections.PathCollection at 0x8eb2650>"
532        ]
533       },
534       {
536        "output_type": "display_data",
538        "text": [
539         "<matplotlib.figure.Figure at 0x864d510>"
540        ]
541       }
542      ],
543      "prompt_number": 125
544     },
545     {
546      "cell_type": "code",
547      "collapsed": false,
548      "input": [
549       "scatter([arctan2(p[0], p[1]) for p in c], [p[0]**2. + p[1]**2. for p in c])\n",
550       "title(\"$r$ vs Angle\")\n",
552       "ylabel(\"$r$\")"
553      ],
554      "language": "python",
556      "outputs": [
557       {
559        "output_type": "pyout",
560        "prompt_number": 126,
561        "text": [
562         "<matplotlib.text.Text at 0x8ea09d0>"
563        ]
564       },
565       {
567        "output_type": "display_data",
569        "text": [
570         "<matplotlib.figure.Figure at 0x8cab9d0>"
571        ]
572       }
573      ],
574      "prompt_number": 126
575     },
576     {
577      "cell_type": "markdown",
579      "source": [
580       "Perfect circle $\\implies r = 1 \\forall \\theta$"
581      ]
582     },
583     {
584      "cell_type": "code",
585      "collapsed": false,
586      "input": [],
587      "language": "python",
589      "outputs": []
590     }
591    ],
593   }
594  ]
595 }

UCC git Repository :: git.ucc.asn.au