Photometric Stereo.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from pybug.io import auto_import\n",
      "from pybug.image import Image\n",
      "import mayavi.mlab as mlab\n",
      "import numpy as np\n",
      "\n",
      "bln = auto_import('/vol/atlas/databases/alex_images/bln*.ppm')\n",
      "# Create a 4 channel image where each channel is the greyscale of an image\n",
      "images = Image(np.concatenate([im.as_greyscale().pixels for im in bln], axis=2))\n",
      "\n",
      "lights = np.array([[ 0.5,  0.4, 2],\n",
      "                   [-0.5,  0.4, 2],\n",
      "                   [-0.5, -0.4, 2],\n",
      "                   [ 0.5, -0.4, 2]])\n",
      "\n",
      "images = images.pixels"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Found 4 files. (0/4) are importable\n",
        "\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Found 4 files. (4/4) are importable\n",
        "\r",
        "Creating importer for <pybug.io.image.PILImporter object at 0x8cc8050> (1 of 4)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        "Creating importer for <pybug.io.image.PILImporter object at 0x8cc8090> (2 of 4)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        "Creating importer for <pybug.io.image.PILImporter object at 0x8cc80d0> (3 of 4)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        "Creating importer for <pybug.io.image.PILImporter object at 0x8cc8110> (4 of 4)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Photometric Stereo"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.linalg import inv, pinv, norm\n",
      "\n",
      "def photometric_stereo(images, lights):\n",
      "    LL = pinv(lights)\n",
      "    \n",
      "    height, width, N = images.shape\n",
      "    albedo = np.zeros([height, width])\n",
      "    \n",
      "    n = np.zeros([height, width, 3])\n",
      "    p = np.zeros([height, width])\n",
      "    q = np.zeros([height, width])\n",
      "    \n",
      "    for ii in xrange(height):\n",
      "        for jj in xrange(width):\n",
      "            I = images[ii, jj, :].flatten()\n",
      "            nn = np.dot(LL, I)\n",
      "            pp = norm(nn)\n",
      "            albedo[ii, jj] = pp\n",
      "            \n",
      "            if pp != 0.0:\n",
      "                # normal = n / albedo \n",
      "                n[ii, jj, :] = nn / pp  \n",
      "                if n[ii, jj, 2] != 0.0:\n",
      "                    p[ii, jj] = n[ii,jj, 0] / n[ii, jj, 2]  # x / z\n",
      "                    q[ii, jj] = n[ii,jj, 1] / n[ii, jj, 2]  # y / z\n",
      "                    \n",
      "    return n, albedo"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n, albedo = photometric_stereo(images, lights)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Frankot-Chellapa integration"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.fftpack import ifftshift, fft2, ifft2\n",
      "\n",
      "def frankotchellappa(dzdx, dzdy):\n",
      "    rows, cols = dzdx.shape\n",
      "     \n",
      "    # The following sets up matrices specifying frequencies in the x and y\n",
      "    # directions corresponding to the Fourier transforms of the gradient\n",
      "    # data.  They range from -0.5 cycles/pixel to + 0.5 cycles/pixel.\n",
      "    # The scaling of this is irrelevant as long as it represents a full\n",
      "    # circle domain. This is functionally equivalent to any constant * pi\n",
      "    row_grid = np.linspace(-0.5, 0.5, rows)\n",
      "    col_grid = np.linspace(-0.5, 0.5, cols)\n",
      "    wx, wy = np.meshgrid(col_grid, row_grid)\n",
      "     \n",
      "    # Quadrant shift to put zero frequency at the appropriate edge\n",
      "    wx = ifftshift(wx)\n",
      "    wy = ifftshift(wy)\n",
      " \n",
      "    # Fourier transforms of gradients\n",
      "    DZDX = fft2(dzdx)\n",
      "    DZDY = fft2(dzdy)\n",
      " \n",
      "    # Integrate in the frequency domain by phase shifting by pi/2 and\n",
      "    # weighting the Fourier coefficients by their frequencies in x and y and\n",
      "    # then dividing by the squared frequency.  eps is added to the\n",
      "    # denominator to avoid division by 0.  \n",
      "    Z = (-1j * wx * DZDX - 1j * wy * DZDY) / (wx ** 2 + wy ** 2 + 0.000001)\n",
      "    \n",
      "    return np.real(ifft2(Z))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 258
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hh = frankotchellappa(n[:, :, 0].copy() / n[:, :, 2].copy(), n[:, :, 1].copy() / n[:, :, 2].copy())\n",
      "mlab.surf(hh, warp_scale='auto')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 259,
       "text": [
        "<mayavi.modules.surface.Surface at 0x7f5adf9a9a70>"
       ]
      }
     ],
     "prompt_number": 259
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Poisson Reconstruction Using Neumann boundary conditions"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def weighted_finite_differences(gx, gy, wx, wy):\n",
      "    # Fix gradients at boundaries to 0 (neumann boundary condition)\n",
      "    gx[:, -1] = 0\n",
      "    gy[-1, :] = 0\n",
      "    \n",
      "    # Weight the gradients\n",
      "    gx = gx * wx\n",
      "    gy = gy * wy\n",
      "    \n",
      "    # Pad each array around both axes by one pixels with the value 0\n",
      "    gx = np.pad(gx, 1, mode='constant')\n",
      "    gy = np.pad(gy, 1, mode='constant')\n",
      "    gxx = np.zeros_like(gx)\n",
      "    gyy = np.zeros_like(gx)\n",
      "    \n",
      "    # Finite differences\n",
      "    # Take the finite differences:\n",
      "    #     gyy[j+1, k] = gy[j+1,k] - gy[j,k]\n",
      "    #     gxx[j, k+1] = gx[j,k+1] - gx[j,k]\n",
      "    # where \n",
      "    #     j = 0:height+1,\n",
      "    #     k = 0:width+1\n",
      "    gyy[1:, :-1] = gy[1:, :-1] - gy[:-1, :-1]\n",
      "    gxx[:-1, 1:] = gx[:-1, 1:] - gx[:-1, :-1]\n",
      "    f = gxx + gyy\n",
      "    f = f[1:-1, 1:-1]\n",
      "    \n",
      "    return f"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.fftpack import dct, idct\n",
      "\n",
      "def poisson_neumann(gx, gy):   \n",
      "    height, width = gx.shape\n",
      "    \n",
      "    f = weighted_finite_differences(gx, gy, np.ones_like(gx), np.ones_like(gy))\n",
      "    \n",
      "    # Compute cosine transform\n",
      "    fcos = dct(dct(f, norm='ortho', axis=0), norm='ortho', axis=1)\n",
      "    \n",
      "    # Compute the solution in the fourier domain\n",
      "    x, y = meshgrid(np.arange(width), np.arange(height))\n",
      "    denom = (2.0 * np.cos(np.pi * x / width) - 2.0) + (2.0 * np.cos(np.pi * y / height) - 2.0)\n",
      "    # First element is fixed to 0.0\n",
      "    fcos /= denom\n",
      "    fcos[0, 0] = 0\n",
      "\n",
      "    # Compute inverse discrete cosine transform to find solution\n",
      "    # in spatial domain\n",
      "    return idct(idct(fcos, norm='ortho', axis=0), norm='ortho', axis=1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hh = poisson_neumann(n[:, :, 0].copy() / n[:, :, 2].copy(), n[:, :, 1].copy() / n[:, :, 2].copy())\n",
      "mlab.surf(hh, warp_scale='auto')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 262,
       "text": [
        "<mayavi.modules.surface.Surface at 0x7f5adf8d17d0>"
       ]
      }
     ],
     "prompt_number": 262
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# M-estimator reconstruction"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.sparse.linalg import spsolve, lsqr\n",
      "from scipy.sparse import csc_matrix, issparse\n",
      "\n",
      "def weighted_laplacian(height, width, wx, wy):\n",
      "    # Laplacian matrix for D*[Zx;Zy] where D is 2*2 diffusion tensor\n",
      "    # kernel D should be symmetric   \n",
      "    \n",
      "    wx = np.pad(wx, 1, mode='constant')\n",
      "    wy = np.pad(wy, 1, mode='constant')\n",
      "    \n",
      "    N = (height + 2) * (width + 2)\n",
      "    mask = np.zeros([height + 2, width + 2], dtype=np.int32)\n",
      "    # Set to all 1s inside the padding\n",
      "    mask[1:-1, 1:-1] = 1\n",
      "    # Equivalent to Matlab's 'find' method (linear indices in to mask matrix)\n",
      "    idx = np.ravel_multi_index((mask == 1).nonzero(), mask.shape)\n",
      "\n",
      "    wx_flat = wx.flatten()\n",
      "    wy_flat = wy.flatten()\n",
      "    \n",
      "    A = csc_matrix((-wx_flat[idx], (idx, idx + height + 2)), shape=(N, N))\n",
      "    A = A + csc_matrix((-wx_flat[idx - height - 2], (idx, idx - height - 2)), shape=(N, N))\n",
      "    A = A + csc_matrix((-wy_flat[idx], (idx, idx + 1)), shape=(N, N))\n",
      "    A = A + csc_matrix((-wy_flat[idx - 1], (idx, idx - 1)), shape=(N, N))\n",
      "    \n",
      "    A = A[np.ix_(idx, idx)]\n",
      "    # width * height\n",
      "    N = A.shape[0]\n",
      "    # Sum over the columns and then turn back in to a 1D array\n",
      "    dd = np.asarray(A.sum(1)).flatten()\n",
      "    idx = np.arange(N)\n",
      "    return A + csc_matrix((-dd, (idx, idx)), shape=(N, N))\n",
      "\n",
      "def huber_weights(k, previous_estimate, xs, ys):\n",
      "    gradient_y, gradient_x = np.gradient(previous_estimate)\n",
      "\n",
      "    # Huber estimate for xs\n",
      "    x_estimate = np.abs(gradient_x - xs)\n",
      "    w_x = np.zeros_like(xs)\n",
      "    idx = x_estimate <= k\n",
      "    w_x[idx] = 1\n",
      "    idx = x_estimate >= k\n",
      "    w_x[idx] = k / x_estimate[idx]\n",
      "    \n",
      "    # Huber estimate for ys\n",
      "    y_estimate = np.abs(gradient_y - ys)\n",
      "    w_y = np.zeros_like(ys)\n",
      "    idy = y_estimate <= k\n",
      "    w_y[idy] = 1\n",
      "    idy = y_estimate >= k\n",
      "    w_y[idy] = k / y_estimate[idy]\n",
      "    \n",
      "    return w_x, w_y\n",
      "\n",
      "def M_estimator(gx, gy, max_iters=3, k=1.345):\n",
      "    height, width = gx.shape\n",
      "    \n",
      "    # Calculate an initial estimate from the poisson solver\n",
      "    previous_estimate = poisson_neumann(gx, gy)\n",
      "    \n",
      "    for ii in xrange(max_iters):\n",
      "        # Huber M-estimator\n",
      "        w_gx, w_gy = huber_weights(k, previous_estimate, gx, gy)\n",
      "    \n",
      "        f = weighted_finite_differences(gx, gy, w_gx, w_gy)\n",
      "        A = weighted_laplacian(height, width, w_gx, w_gy)\n",
      "    \n",
      "        current_estimate = spsolve(-A, f.flatten())\n",
      "#         current_estimate = np.concatenate([[0.0], current_estimate])\n",
      "        current_estimate = np.reshape(current_estimate, [height, width])\n",
      "    \n",
      "        idx = np.abs(previous_estimate) > 0.01\n",
      "        if(idx.size):\n",
      "            ee = (previous_estimate[idx] - current_estimate[idx]) / previous_estimate[idx]\n",
      "            ee = ee ** 2\n",
      "            ee = np.mean(ee)\n",
      "            if ee < 1e-6:\n",
      "                break\n",
      "                \n",
      "        previous_estimate = current_estimate\n",
      "    \n",
      "    return current_estimate"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 32
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hh = M_estimator(n[:, :, 0].copy() / n[:, :, 2].copy(), n[:, :, 1].copy() / n[:, :, 2].copy())\n",
      "mlab.surf(hh, warp_scale='auto')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[-204.77693239 -204.22408351 -203.62130444 ..., -176.69001396\n",
        "  -176.98499319 -177.32890847]\n",
        " [-180.08791029 -180.4666254  -180.93251607 ..., -164.2616689  -163.73226038\n",
        "  -163.24163583]\n",
        " [-159.89439392 -159.24225381 -158.51932677 ..., -184.93173493\n",
        "  -185.55603577 -186.23112722]\n",
        " ..., \n",
        " [  40.19136702   41.15371148   42.16608469 ...,  127.79005948\n",
        "   129.21519604  130.16336329]\n",
        " [ 134.99619513  135.36448789  135.61856195 ...,  128.10785393\n",
        "   129.05372199  129.45948195]\n",
        " [ 131.23634806  131.38425187  131.59593912 ...,  175.05628793\n",
        "   174.61382538  174.23159411]]\n",
        "[[-1095.01222534 -1083.03955235 -1071.34542216 ...,  1635.19353371\n",
        "   1607.62540209  1579.56718492]\n",
        " [ 1469.12324145  1441.86897822  1412.39170695 ..., -1109.48057189\n",
        "  -1098.90305286 -1088.26290198]\n",
        " [-1038.65596884 -1027.03575168 -1014.12212247 ...,  1893.43331163\n",
        "   1861.01210233  1826.82590008]\n",
        " ..., \n",
        " [ -555.58095963  -587.17352524  -620.64620286 ...,  -274.83722157\n",
        "   -250.0706466   -228.14177174]\n",
        " [ -112.27846078   -82.622408     -55.65385259 ...,   -74.30806973\n",
        "    -85.04805129   -98.63077658]\n",
        " [ -189.21314262  -206.25181422  -225.70473595 ...,   888.27791772\n",
        "    922.63593896   952.17619055]]"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "[[-2613.08002722 -2701.98562589 -2792.45318914 ...,  -121.2189113\n",
        "    -19.53972412    81.30307136]\n",
        " [  456.3934514    510.00559493   581.43079992 ..., -2293.0605572\n",
        "  -2384.60097755 -2470.80151339]\n",
        " [-2793.38733824 -2868.44087992 -2948.33209423 ..., -3435.78618655\n",
        "  -3335.10217228 -3254.61759266]\n",
        " ..., \n",
        " [ 3241.57305324  3260.55972923  3272.1382633  ..., -2599.12794785\n",
        "  -2600.6857526  -2603.28095768]\n",
        " [-2612.48992592 -2618.37113924 -2624.19398782 ...,  4111.9589537\n",
        "   4152.2868991   4186.19682657]\n",
        " [ 4268.79772605  4284.75227262  4309.45237056 ..., -3141.55098604\n",
        "  -3148.35179488 -3156.61226294]]"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 33,
       "text": [
        "<mayavi.modules.surface.Surface at 0x10381050>"
       ]
      }
     ],
     "prompt_number": 33
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.io import loadmat\n",
      "\n",
      "H = 10\n",
      "W = 10\n",
      "gx = np.ones([10, 10])\n",
      "gy = np.ones([10, 10])\n",
      "k = 1.345\n",
      "wx = loadmat('/home/pts08/Desktop/matlab.mat')['wx']\n",
      "wy = loadmat('/home/pts08/Desktop/matlab.mat')['wy']\n",
      "\n",
      "g_wx, g_wy = huber_weights(k, wx, gx, gy)\n",
      "print g_wx\n",
      "print g_wy"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 0.8117  1.      1.      1.      1.      1.      1.      1.      0.9729  1.    ]\n",
        " [ 1.      0.9373  0.9154  1.      1.      1.      1.      1.      1.      0.7569]\n",
        " [ 1.      1.      1.      1.      1.      1.      1.      1.      1.      1.    ]\n",
        " [ 0.9419  1.      1.      1.      1.      1.      1.      1.      1.      1.    ]\n",
        " [ 1.      1.      0.9951  1.      1.      1.      1.      1.      1.      1.    ]\n",
        " [ 1.      1.      1.      1.      1.      1.      1.      1.      1.      1.    ]\n",
        " [ 1.      1.      1.      1.      1.      1.      0.9735  1.      1.      1.    ]\n",
        " [ 1.      1.      1.      1.      1.      1.      1.      1.      1.      1.    ]\n",
        " [ 1.      1.      1.      1.      1.      1.      1.      1.      1.      0.9828]\n",
        " [ 1.      0.9629  0.9196  1.      1.      1.      1.      1.      1.      1.    ]]\n",
        "[[ 1.      1.      0.8302  0.8034  1.      1.      0.899   0.8478  1.      1.    ]\n",
        " [ 1.      1.      1.      1.      1.      1.      1.      1.      1.      1.    ]\n",
        " [ 1.      1.      1.      1.      1.      1.      1.      1.      1.      1.    ]\n",
        " [ 1.      1.      1.      1.      1.      1.      1.      1.      1.      1.    ]\n",
        " [ 0.9553  1.      1.      1.      1.      1.      1.      1.      1.      1.    ]\n",
        " [ 1.      1.      1.      1.      1.      1.      1.      0.9844  1.      1.    ]\n",
        " [ 1.      1.      1.      1.      1.      1.      0.9537  1.      1.      1.    ]\n",
        " [ 1.      1.      1.      1.      1.      1.      1.      1.      1.      1.    ]\n",
        " [ 1.      1.      1.      1.      1.      1.      1.      1.      1.      1.    ]\n",
        " [ 1.      1.      0.9062  0.7021  1.      0.9879  1.      1.      1.      1.    ]]\n"
       ]
      }
     ],
     "prompt_number": 27
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy.linalg import eig\n",
      "from scipy.ndimage.filters import gaussian_filter\n",
      "\n",
      "def diffusion_kernel_function(gx, gy):\n",
      "\n",
      "    # Tensor based\n",
      "    # Instead of individual weights to gx and gy, we have a 2*2 tensor getting\n",
      "    # multiplied by [gx,gy] vector at each pixel\n",
      "    # get tensor values using local analysis\n",
      "    \n",
      "    H, W = gx.shape\n",
      "    p = gx\n",
      "    q = gy\n",
      "    \n",
      "    sigma = 0.5\n",
      "# How do we create a kernel size in Scipy?\n",
      "#     ss = np.floor(6.0 * sigma)\n",
      "#     if(ss < 3):\n",
      "#         ss = 3\n",
      "#     ww = fspecial('gaussian',ss,sigma)\n",
      "    \n",
      "    \n",
      "    T11 = gaussian_filter(p ** 2, sigma)\n",
      "    T22 = gaussian_filter(q ** 2, sigma)\n",
      "    T12 = gaussian_filter(p * q, sigma)\n",
      "    \n",
      "    # find eigen values\n",
      "    ImagPart = np.sqrt((T11 - T22) ** 2 + 4 * (T12 ** 2));\n",
      "    EigD_1 = (T22 + T11 + ImagPart) / 2.0\n",
      "    EigD_2 = (T22 + T11 - ImagPart) / 2.0\n",
      "    \n",
      "    alpha = 0.02\n",
      "    THRESHOLD_SMALL = np.max(EigD_1.flatten()) / 100.0\n",
      "    \n",
      "    L1 = np.ones([H, W]);\n",
      "    idx = EigD_1 > THRESHOLD_SMALL\n",
      "    L1[idx] = alpha + 1.0 - np.exp(-3.315 / (EigD_1[idx] ** 4))\n",
      "    L2 = np.ones([H, W])\n",
      "    \n",
      "    D11 = np.zeros([H, W])\n",
      "    D12 = np.zeros([H, W])\n",
      "    D22 = np.zeros([H, W])\n",
      "\n",
      "    for ii in xrange(H):\n",
      "        for jj in xrange(W):\n",
      "            Wmat = np.array([[T11[ii, jj], T12[ii, jj]], \n",
      "                             [T12[ii, jj], T22[ii, jj]]])\n",
      "            v, d = eig(Wmat)\n",
      "            v = np.diagflat(v)\n",
      "            if d[0, 0] > d[1, 1]:\n",
      "                d0 = d[0, 0]\n",
      "                d[0, 0] = d[1, 1]\n",
      "                d[1, 1] = d0\n",
      "                \n",
      "                v0 = v1[:, 0].copy()\n",
      "                v1[:, 0] = v[:, 1]\n",
      "                v1[:, 1] = v0\n",
      "    \n",
      "            # d(1,1) is smaller\n",
      "            d[0, 0] = L2[ii, jj]\n",
      "            d[1, 1] = L1[ii, jj]\n",
      "    \n",
      "            Wmat = np.dot(np.dot(v, d), v.T)\n",
      "            D11[ii, jj] = Wmat[0, 0]\n",
      "            D22[ii, jj] = Wmat[1, 1]\n",
      "            D12[ii, jj] = Wmat[0, 1]\n",
      "    \n",
      "    A = laplacian_matrix_tensor(H,W,D11,D12,D12,D22)\n",
      "    f = calculate_f_tensor(p,q,D11,D12,D12,D22)\n",
      "    \n",
      "    Z = spsolve(A, f.flatten())\n",
      "#     Z = [0;Z];\n",
      "    Z = np.reshape(Z, [H,W])\n",
      "    Z = Z - np.min(Z)\n",
      "    return Z"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hh = diffusion_kernel_function(n[:, :, 0].copy(), n[:, :, 1].copy())\n",
      "mlab.surf(hh, warp_scale='auto')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "<mayavi.modules.surface.Surface at 0xedc7f50>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}