author Sam Moore Thu, 16 Jan 2014 05:36:01 +0000 (13:36 +0800) committer Sam Moore Thu, 16 Jan 2014 05:36:01 +0000 (13:36 +0800)
Using Goldman's Iterated Function System to generate a Bezier is horribly inefficient.
I guess the point is that it proves that a Bezier is a fractal.

The subdivision still looks like it gives a slightly different curve to
just sampling P(t), but maybe I aren't using enough points.

index 976a6b6..c7a2805 100644 (file)
"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",
+      "I think the algorithm developed by Goldman might be worse than just applying the Bezier definition...\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"
+      "Binomial Coefficients"
]
},
{
@@ -78,7 +56,7 @@
"language": "python",
"outputs": [],
-     "prompt_number": 216
+     "prompt_number": 2
},
{
"cell_type": "code",
"output_type": "stream",
"stream": "stdout",
"text": [
-        "0.000294923782349\n",
-        "0.000556945800781\n"
+        "0.00108003616333\n",
+        "0.000107049942017\n"
]
}
],
-     "prompt_number": 215
+     "prompt_number": 3
},
{
"Bezier from 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)!} = \\text{nCj}$"
+     ]
+    },
{
"cell_type": "code",
"collapsed": false,
"    \"\"\" 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 BezierCurve(P, nPoints):\n",
+      "    return map(lambda t : Bezier(P,t), linspace(0,1,nPoints))\n",
"        \n",
-      "def PlotBezier(P,nPoints=50,style='o-'):\n",
+      "def PlotBezier(P,nPoints=50,style='-'):\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",
"language": "python",
"outputs": [],
-     "prompt_number": 266
+     "prompt_number": 4
},
{
"cell_type": "code",
"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()"
+      "PlotBezier([(0,0),(0.4,-0.8),(1,1)],nPoints=100,style='-')"
],
"language": "python",
"outputs": [
+      {
+       "output_type": "stream",
+       "stream": "stderr",
+       "text": [
+        "/usr/lib/pymodules/python2.7/matplotlib/axes.py:4486: UserWarning: No labeled objects found. Use label='...' kwarg on individual plots.\n",
+        "  warnings.warn(\"No labeled objects found. \"\n"
+       ]
+      },
{
"output_type": "display_data",
"text": [
-        "<matplotlib.figure.Figure at 0x3c4ee50>"
+        "<matplotlib.figure.Figure at 0x2bb2bd0>"
]
}
],
-     "prompt_number": 267
+     "prompt_number": 5
},
{
"level": 2,
"source": [
-      "Fractal Beziers"
+      "De Casteljau Algorithm"
]
},
{
"cell_type": "markdown",
"source": [
-      "$L P = Q$ and $M P = R$"
+      "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",
+      "\n",
+      "Defining $(n+1)\\times(n+1)$ matrices $L$ and $M$ such that $L P = Q$ and $M P = R$\n",
+      "\n",
+      "Then it can be shown (subdividing at $t = 1/2$)\n",
+      "\n",
+      "$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)$"
+     ]
+    },
+    {
+     "cell_type": "markdown",
+     "source": [
+      "**Confusion:** Do $Q$ and $R$ converge to the same curve as applying the definition to produce points along $P(t)$?\n",
+      "\n",
+      "It is implied in Goldman, but it doesn't look like it, unless I made a mistake. See below."
]
},
{
"collapsed": false,
"input": [
"def L(n):\n",
+      "    \"\"\" Gives L of size (n+1)*(n+1) \"\"\"\n",
"    try:\n",
"        return L.dynamic[n]\n",
"    except KeyError:\n",
"L.dynamic = {}                        \n",
"\n",
"def M(n):\n",
+      "    \"\"\" Gives M of size (n+1)*(n+1) \"\"\"\n",
"    try:\n",
"        return M.dynamic[n]\n",
"    except KeyError:\n",
"language": "python",
"outputs": [],
-     "prompt_number": 296
+     "prompt_number": 6
},
{
"cell_type": "code",
"collapsed": false,
"input": [
-      "points = [(0,0), (0.5,1), (1,0)]"
+      "L(2)"
],
"language": "python",
-     "outputs": [],
-     "prompt_number": 298
+     "outputs": [
+      {
+       "output_type": "pyout",
+       "prompt_number": 7,
+       "text": [
+        "array([[ 1.  ,  0.  ,  0.  ],\n",
+        "       [ 0.5 ,  0.5 ,  0.  ],\n",
+        "       [ 0.25,  0.5 ,  0.25]])"
+       ]
+      }
+     ],
+     "prompt_number": 7
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def Casteljau(P, nPoints):\n",
+      "    \"\"\" Approximates the Bezier curve with control points P using the de Casteljau subdivision algorithm \"\"\"\n",
"    while len(P) < nPoints:\n",
"        Q = dot(L(len(P)-1), P)\n",
"        R = dot(M(len(P)-1), P)\n",
-      "        P = list(Q) + list(R)[1:]\n",
+      "        P = list(Q) + list(R)\n",
"    return P\n",
"    "
],
"language": "python",
"outputs": [],
-     "prompt_number": 329
+     "prompt_number": 8
},
{
"cell_type": "code",
"collapsed": false,
"input": [
-      "points = [(0,0), (0.5,1), (1,0)]\n",
+      "points = [(0,0), (0.5,1),(1,0),(1.4,-2)]\n",
"c = Casteljau(points,50)\n",
"PlotBezier(points)\n",
-      "plot([p[0] for p in c], [p[1] for p in c], 'go-')"
+      "plot([p[0] for p in c], [p[1] for p in c], 'g-')"
],
"language": "python",
{
"output_type": "pyout",
-       "prompt_number": 332,
+       "prompt_number": 9,
"text": [
-        "[<matplotlib.lines.Line2D at 0x5a3a6d0>]"
+        "[<matplotlib.lines.Line2D at 0x2af9b10>]"
]
},
{
"output_type": "display_data",
"text": [
-        "<matplotlib.figure.Figure at 0x5942790>"
+        "<matplotlib.figure.Figure at 0x2af9950>"
]
}
],
-     "prompt_number": 332
+     "prompt_number": 9
},
{
"cell_type": "markdown",
"source": [
-      "So, it sort of works but not giving the same curve as the vanilla method"
+      "So, it sort of works but not giving the same curve? Or maybe I just need *lots* more points and they will converge."
+     ]
+    },
+    {
+     "level": 2,
+     "source": [
+      "Iterated Function Systems"
+     ]
+    },
+    {
+     "cell_type": "markdown",
+     "source": [
+      "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",
+      "\n",
+      "\n",
+      "$\\{L_p, M_p\\}$ is an iterated function system, $\\implies$ Beziers are fractals.\n",
+      "\n",
+      "Problem: $P$ needs to be invertable; introduce homogeneous coordinates $\\langle$insert hand waving here$\\rangle$ to \"lift\" it and then it all magically works"
]
},
{
"cell_type": "code",
"collapsed": false,
-     "input": [],
+     "input": [
+      "def lift(P):\n",
+      "    \"\"\" Lift P to make it square by introducing homogeneous coordinates\"\"\"\n",
+      "    n = len(P)\n",
+      "    result = []\n",
+      "    for i in xrange(len(P)):\n",
+      "        ones = min(n-len(P[i]), max(1, i-1))\n",
+      "        zeroes = max(n-len(P[i])-ones, 0)\n",
+      "        result += [asarray(list(P[i]) + [1. for _ in xrange(ones)] + [0. for _ in xrange(zeroes)])]\n",
+      "    return asarray(result)\n",
+      "\n",
+      "def lower(P):\n",
+      "    \"\"\" Lower square P back to list of 2D points \"\"\"\n",
+      "    return [row[0:2] for row in P]\n",
+      "\n",
+      "def Lp(P):\n",
+      "    return dot(inv(P), dot(L(len(P)-1), P))\n",
+      "\n",
+      "def Mp(P):\n",
+      "    return dot(inv(P), dot(M(len(P)-1), P)) # Typo in Goldman? Reads \"P^-1 * R * P\" instead of \"P^-1 * M * P\""
+     ],
"language": "python",
-     "outputs": []
+     "outputs": [],
+     "prompt_number": 10
+    },
+    {
+     "cell_type": "markdown",
+     "source": [
+      "So, Goldman's method is:\n",
+      "\n",
+      "<ol>\n",
+      " <li> lift the control points\n",
+      "P from the ambient affine space of dimension $d$ to the\n",
+      "ambient affine space of dimension $n$;\n",
+      " </li><li> generate the corresponding higher dimensional Bezier curve using the\n",
+      "iterated function system $\\{L_p , M_p \\}$ ; \n",
+      "</li><li> project the\n",
+      "resulting $n$-dimensional Bezier curve orthogonally back\n",
+      "down to the original dimension $d$.\n",
+      "</ol>\n"
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "def Goldman(P, nPoints):\n",
+      "    while len(P) < nPoints:\n",
+      "        Pl = lift(P)\n",
+      "        LP = Lp(Pl)\n",
+      "        MP = Mp(Pl)\n",
+      "        Q = dot(Pl, LP)\n",
+      "        R = dot(Pl, MP)\n",
+      "        P = lower(Q) + lower(R)\n",
+      "    return P\n",
+      "        "
+     ],
+     "language": "python",
+     "outputs": [],
+     "prompt_number": 11
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "points = [(0,0), (0.5,1),(1,0),(1.4,-2)]\n",
+      "c = Casteljau(points,50)\n",
+      "c2 = Goldman(points,50)\n",
+      "PlotBezier(points)\n",
+      "plot([p[0] for p in c], [p[1] for p in c], 'g-')\n",
+      "plot([p[0] for p in c2], [p[1] for p in c2], 'r-')"
+     ],
+     "language": "python",
+     "outputs": [
+      {
+       "output_type": "pyout",
+       "prompt_number": 12,
+       "text": [
+        "[<matplotlib.lines.Line2D at 0x30a8750>]"
+       ]
+      },
+      {
+       "output_type": "display_data",
+       "text": [
+        "<matplotlib.figure.Figure at 0x2bfe7d0>"
+       ]
+      }
+     ],
+     "prompt_number": 12
+    },
+    {
+     "cell_type": "markdown",
+     "source": [
+      "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."
+     ]
+    },
+    {
+     "level": 2,
+     "source": [
+      "Time Performance of Casteljau, Goldman, and Basic Bezier creation"
+     ]
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "from time import time\n",
+      "points = [(0,0), (0.5,1),(1,0),(1.4,-2)]\n",
+      "costs = {BezierCurve : [], Casteljau : [], Goldman : []}\n",
+      "stddev = {}\n",
+      "nAverages = 10\n",
+      "depth=8\n",
+      "for method in costs:\n",
+      "\n",
+      "    stddev.update({method : []})\n",
+      "    for n in xrange(1,depth):\n",
+      "        #display(\"depth %d for %s\" % (n, str(method)))\n",
+      "        c = []\n",
+      "        for _ in xrange(nAverages):\n",
+      "            t0 = time()\n",
+      "            p = method(points, nPoints=len(points) * 2.**n)\n",
+      "            c += [time()-t0]\n",
+      "        costs[method] += [mean(c)]\n",
+      "        stddev[method] += [var(c)**0.5]"
+     ],
+     "language": "python",
+     "outputs": [],
+     "prompt_number": 17
+    },
+    {
+     "cell_type": "code",
+     "collapsed": false,
+     "input": [
+      "title(\"Time to produce $n$-dimensional Bezier Curve\")\n",
+      "xlabel(\"$n$\")\n",
+      "ylabel(\"$t$ (s)\")\n",
+      "x = [len(points) * 2.**n for n in xrange(1,depth)]\n",
+      "for method in costs:\n",
+      "    errorbar(x,costs[method],yerr=stddev[method],label =str(method).split(\" \")[1])\n",
+      "legend(loc=\"upper left\")"
+     ],
+     "language": "python",
+     "outputs": [
+      {
+       "output_type": "pyout",
+       "prompt_number": 18,
+       "text": [
+        "<matplotlib.legend.Legend at 0x2f7f790>"
+       ]
+      },
+      {
+       "output_type": "display_data",
+       "text": [
+        "<matplotlib.figure.Figure at 0x3a1be50>"
+       ]
+      }
+     ],
+     "prompt_number": 18
}
],