author Sam Moore Wed, 15 Jan 2014 14:36:46 +0000 (22:36 +0800) committer Sam Moore Wed, 15 Jan 2014 14:36:46 +0000 (22:36 +0800)
For your information, Potter, a Bezier is a stone taken from the stomache of a goat and it will save you from most -

Wait, it's actually a fractal, nevermind.

(The way I made it isn't fractally yet. That's the de Casteljau Algorithm. Will try that next)

 figures/koch_snowflake_n=5.png [new file with mode: 0644] patch | blob figures/stierpinski_tri_n=10000.png [new file with mode: 0644] patch | blob ipython_notebooks/de_Casteljau.ipynb [new file with mode: 0644] patch | blob ipython_notebooks/fractals_basic.ipynb patch | blob | history

diff --git a/figures/koch_snowflake_n=5.png b/figures/koch_snowflake_n=5.png
new file mode 100644 (file)
index 0000000..38785f1
Binary files /dev/null and b/figures/koch_snowflake_n=5.png differ
diff --git a/figures/stierpinski_tri_n=10000.png b/figures/stierpinski_tri_n=10000.png
new file mode 100644 (file)
index 0000000..a8b7986
Binary files /dev/null and b/figures/stierpinski_tri_n=10000.png differ
diff --git a/ipython_notebooks/de_Casteljau.ipynb b/ipython_notebooks/de_Casteljau.ipynb
new file mode 100644 (file)
index 0000000..7e3dbd3
--- /dev/null
@@ -0,0 +1,194 @@
+{
+  "name": ""
+ },
+ "nbformat": 3,
+ "nbformat_minor": 0,
+ "worksheets": [
+  {
+   "cells": [
+    {
+     "cell_type": "markdown",
+     "source": [
+      "Described in Goldman, used as argument that Beziers are fractals because they are fixed points of an iterated function system.\n",
+      "\n",
+      "The de Casteljau algorithm splits a Bezier curve into two Bezier segments.\n",
+      "\n"
+     ]
+    },
+    {
+     "level": 3,
+     "source": [
+      "Bezier Curve Definition"
+     ]
+    },
+    {
+     "cell_type": "markdown",
+     "source": [
+      "\n",
+      "Bezier curve $P(t)$ of degree $n$\n",
+      "\n",
+      "$P(t) = \\sum_{j=0}^{n} B_j^n (t) P_j \\quad \\quad \\quad 0 \\leq t \\leq 1$\n",
+      "\n",
+      "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",
+      "\n",
+      "Binomial Coefficients: $\\left(^n_j\\right) = \\frac{n!}{n!(n-j)!}$"
+     ]
+    },
+    {
+     "level": 2,
+     "source": [
+      "Helpers"
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "from operator import mul    # or mul=lambda x,y:x*y\n",
+      "from fractions import Fraction\n",
+      "\n",
+      "def nCk(n,k): \n",
+      "    \"\"\" Binomial Coefficient nCk \"\"\"\n",
+      "    # Based on:\n",
+      "    # http://stackoverflow.com/questions/3025162/statistics-combinations-in-python\n",
+      "    # http://stackoverflow.com/questions/279561/what-is-the-python-equivalent-of-static-variables-inside-a-function\n",
+      "\n",
+      "    \n",
+      "    # Dynamic programming is possible in python!\n",
+      "    # Not sure if it's actually worth it though\n",
+      "    try:\n",
+      "        return nCk.dynamic[(n,k)]\n",
+      "    except KeyError:\n",
+      "        nCk.dynamic.update({(n,k) : int( reduce(mul, (Fraction(n-i, i+1) for i in range(k)), 1) )})\n",
+      "    return nCk.dynamic[(n,k)]\n",
+      "nCk.dynamic = {}\n",
+      "\n",
+      "# Simpler version...\n",
+      "def nCr(n,r):\n",
+      "    f = math.factorial\n",
+      "    return f(n) / f(r) / f(n-r)"
+     ],
+     "language": "python",
+     "outputs": [],
+     "prompt_number": 216
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "from time import time\n",
+      "t0 = time()\n",
+      "nCk(100,50)\n",
+      "print str(time() - t0)\n",
+      "t0 = time()\n",
+      "nCr(100,50)\n",
+      "print str(time() - t0)"
+     ],
+     "language": "python",
+     "outputs": [
+      {
+       "output_type": "stream",
+       "stream": "stdout",
+       "text": [
+        "0.000294923782349\n",
+        "0.000556945800781\n"
+       ]
+      }
+     ],
+     "prompt_number": 215
+    },
+    {
+     "level": 2,
+     "source": [
+      "Bezier from Definition"
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "def B(n,j, t):\n",
+      "    \"\"\" Bernstein Basis Function of degree n \"\"\"\n",
+      "    assert(0 <= j <= n)\n",
+      "    return nCk(n,j) * t**j * (1. - t)**(n-j)\n",
+      "\n",
+      "def Bezier(P, t):\n",
+      "    \"\"\" Apply Bezier definition to produce a point on the Bezier curve with control points P \"\"\"\n",
+      "    n = len(P)-1\n",
+      "    return [sum([B(n, j, t) * P[j][i] for j in xrange(n+1)]) for i in xrange(len(P[0]))]\n",
+      "        \n",
+      "def PlotBezier(P,nPoints=50,style='o-'):\n",
+      "    \"\"\" Plot a Bezier \"\"\"\n",
+      "    points = map(lambda t : Bezier(P,t), linspace(0,1,nPoints))\n",
+      "    x = [p[0] for p in points]\n",
+      "    y = [p[1] for p in points]\n",
+      "    plot(x,y, style)\n",
+      "    plot([p[0] for p in P], [p[1] for p in P],'ro--')\n"
+     ],
+     "language": "python",
+     "outputs": [],
+     "prompt_number": 266
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "title(\"Bezier Curves of Degree 2\")\n",
+      "PlotBezier([(0,0),(0.4,0.8),(1,1)],nPoints=100,style='-')\n",
+      "PlotBezier([(0,0),(0.4,-0.8),(1,1)],nPoints=100,style='-')\n",
+      "legend()"
+     ],
+     "language": "python",
+     "outputs": [
+      {
+       "output_type": "display_data",
+       "text": [
+        "<matplotlib.figure.Figure at 0x3c4ee50>"
+       ]
+      }
+     ],
+     "prompt_number": 267
+    },
+    {
+     "level": 2,
+     "source": [
+      "Fractal Beziers"
+     ]
+    },
+    {
+     "cell_type": "markdown",
+     "source": [
+      "Urgh"
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [],
+     "language": "python",
+     "outputs": []
+    }
+   ],
+  }
+ ]
+}
\ No newline at end of file
@@ -88,7 +88,7 @@
"language": "python",
"outputs": [],
-     "prompt_number": 193
+     "prompt_number": 436
},
{
"cell_type": "code",
@@ -99,7 +99,7 @@
"language": "python",
"outputs": [],
-     "prompt_number": 194
+     "prompt_number": 437
},
{
"cell_type": "code",
"language": "python",
"outputs": [],
-     "prompt_number": 195
+     "prompt_number": 438
},
{
"cell_type": "code",
"collapsed": false,
"input": [
+      "title(\"Stierpinski Triangles (10000 points)\")\n",
"scatter(x, y, marker=',', s = 1e-3)"
],
"language": "python",
{
"output_type": "pyout",
-       "prompt_number": 196,
+       "prompt_number": 440,
"text": [
-        "<matplotlib.collections.PathCollection at 0x79a6290>"
+        "<matplotlib.collections.PathCollection at 0x11957d50>"
]
},
{
"output_type": "display_data",
"text": [
-        "<matplotlib.figure.Figure at 0x7565150>"
+        "<matplotlib.figure.Figure at 0x110c8590>"
]
}
],
-     "prompt_number": 196
+     "prompt_number": 440
},
{
"    normal = asarray([-s[1], s[0]])\n",
"    normal = normal / sum(normal*normal)**0.5\n",
"    # Calculate the three points \n",
-      "    # (I *think* that in the definition these are the maps that form the iterated function system)\n",
+      "    # Are these are the maps that form the iterated function system?\n",
+      "    # \"Notice that in each case the entire fractal can be recovered from a single point.\" Goldman\n",
+      "    # But for this I need *two* points :S\n",
+      "    \n",
"    a = asarray(x0) + s/3. # 1/3 from x0 to x1\n",
"    b = asarray(x0) + 2.*s/3. # 2/3 from x0 to x1\n",
"    c = asarray(x0) + 0.5*s + ((sum(s*s)**0.5)/3.)*normal # Form an equilateral triangle with a & b\n",
"    points = koch(n, [(0,0),(0.5, 0.75), (1,0)])\n",
"    x = [p[0] for p in points]\n",
"    y = [p[1] for p in points]\n",
-      "    plot(x, y, ',-')"
+      "    plot(x, y, ',-')\n",
+      "title(\"Koch Snowflake (first 5 iterations)\")"
],
"language": "python",
"outputs": [
+      {
+       "output_type": "pyout",
+       "prompt_number": 434,
+       "text": [
+        "<matplotlib.text.Text at 0x17e28a50>"
+       ]
+      },
{
"output_type": "display_data",
"text": [
-        "<matplotlib.figure.Figure at 0xdd3a110>"
+        "<matplotlib.figure.Figure at 0x110e3310>"
]
}
],
-     "prompt_number": 373
+     "prompt_number": 434
},
{
"cell_type": "code",
"cell_type": "markdown",
"source": [
-      "So I worked it out, but I didn't work it out very efficiently.\n",
-      "It appears to be $O(n^4)$ :S"
+      "So I worked it out, but I didn't work it out very efficiently..."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
+      "# It might be O(n^4) if I remember DSA right which I might not...\n",
"plot(range(nMax), costs)\n",
-      "for deg in xrange(4,5):\n",
-      "    fit = numpy.polynomial.polynomial.polyfit(range(nMax), costs, deg)\n",
-      "    plot(linspace(0, nMax), map(lambda x : sum([fit[i]*x**i for i in xrange(len(fit))]), linspace(0, nMax)))"
+      "fit = numpy.polynomial.polynomial.polyfit(range(nMax), costs, 4)\n",
+      "plot(linspace(0, nMax), map(lambda x : sum([fit[i]*x**i for i in xrange(len(fit))]), linspace(0, nMax)))"
],
"language": "python",
"outputs": [
+      {
+       "output_type": "pyout",
+       "prompt_number": 447,
+       "text": [
+        "[<matplotlib.lines.Line2D at 0x12142190>]"
+       ]
+      },
{
"output_type": "display_data",