rankine-delaunay.ipynb

{
 "metadata": {
  "name": "rankine-delaunay"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "\n",
      "from __future__ import division\n",
      "import numpy as np\n",
      "import matplotlib as mpl\n",
      "import matplotlib.pyplot as plt\n",
      "from scipy.spatial import Delaunay\n",
      "\n",
      "def triangulate_points(points):\n",
      "    \"\"\"Triangulate a set of data.  Possible formats include 'list', 'pairs' and ''.\"\"\"\n",
      "    points = np.array(points)\n",
      "    triangles = Delaunay(points)\n",
      "    return triangles\n",
      "\n",
      "def plot_triangulated_points(points,triangles):\n",
      "    points = np.array(points)\n",
      "    plt.triplot(points[:,0], points[:,1], triangles.vertices.copy(), 'k-')\n",
      "    plt.plot(points[:,0], points[:,1], 'b.')\n",
      "    plt.show()\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "ERROR: Line magic function `%matplotlib` not found."
       ]
      }
     ],
     "prompt_number": 94
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def dot(pt1, pt2):\n",
      "    return (pt1[0]*pt2[0] + pt1[1]*pt2[1])\n",
      "\n",
      "def calc_bary_coords(pt, pta, ptb, ptc):\n",
      "    \"\"\"Calculates the barycentric coordinates of the pair pt in the triangle of pairs pta, ptb, and ptc.\"\"\"\n",
      "    # Move the origin to pta.\n",
      "    vc = [None, None]\n",
      "    vc[0] = ptc[0] - pta[0]; vc[1] = ptc[1] - pta[1] #vc is from pta to ptc\n",
      "    vb = [None, None]\n",
      "    vb[0] = ptb[0] - pta[0]; vb[1] = ptb[1] - pta[1] #vb is from pta to ptb\n",
      "    vp = [None, None]\n",
      "    vp[0] = pt[0] - pta[0];  vp[1] = pt[1] - pta[1]  #vt is from pta to pt\n",
      "    \n",
      "    # Compute the dot products.\n",
      "    dot_cc = dot(vc, vc)\n",
      "    dot_cb = dot(vc, vb)\n",
      "    dot_cp = dot(vc, vp)\n",
      "    dot_bb = dot(vb, vb)\n",
      "    dot_bp = dot(vb, vp)\n",
      "    \n",
      "    # Compute barycentric coordinates\n",
      "    invDenom = 1 / (dot_cc * dot_bb - dot_cb * dot_cb)\n",
      "    u = (dot_bb * dot_cp - dot_cb * dot_bp) * invDenom\n",
      "    v = (dot_cc * dot_bp - dot_cb * dot_cp) * invDenom\n",
      "    \n",
      "    # Check if point is in triangle.\n",
      "    in_tri = ((u >= 0) and (v >= 0) and (u + v < 1))\n",
      "    return ([u, v], in_tri)\n",
      " "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 137
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Retrieve temperature data from online source.\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Generate random data points.\n",
      "from random import random\n",
      "hi = 5\n",
      "sites = [] # sites are the x-y coordinates of the observed temperatures\n",
      "for i in range(0,hi):\n",
      "    site = [random()]\n",
      "    site.append(random())\n",
      "    sites.append(site)\n",
      "\n",
      "triangles = triangulate_points(sites)\n",
      "plot_triangulated_points(sites,triangles)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 156
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pts = [[0.5,0.5]] # pts are the x-y coordinates of the desired site(s)\n",
      "pts_tri = [None]*len(pts)\n",
      "\n",
      "for i, pt in enumerate(pts):\n",
      "    # For a given interpolant point, find the triangle in which it lies.\n",
      "    for j, triangle in enumerate(triangles.vertices):\n",
      "        a = points[triangle[0]]\n",
      "        b = points[triangle[1]]\n",
      "        c = points[triangle[2]]\n",
      "        \n",
      "        # Calculate the barycentric coordinates of an interior point.\n",
      "        (pt_bc,in_tri) = calc_bary_coords(pt,a,b,c)\n",
      "        \n",
      "        if in_tri:\n",
      "            pts_tri[i] = j\n",
      "            break\n",
      "\n",
      "#for i in range(0,len(pts)):\n",
      "#    print(pts[i], pts_tri[i], triangles.vertices[i],points[triangles.vertices[pts_tri[i]][0]],points[triangles.vertices[pts_tri[i]][1]],points[triangles.vertices[pts_tri[i]][2]])\n",
      "\n",
      "sites = np.array(sites)\n",
      "pts = np.array(pts)\n",
      "plt.triplot(sites[:,0], sites[:,1], triangles.vertices.copy(), 'k-')\n",
      "plt.plot(sites[:,0], sites[:,1], 'b.')\n",
      "plt.plot(pts[:,0], pts[:,1], 'g.')\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 167
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(triangles.vertices)\n",
      "print(sites)\n",
      "\n",
      "for i in range(0,len(pts)):\n",
      "    print(pts[i], pts_tri[i], triangles.vertices[pts_tri[i]], sites[triangles.vertices[pts_tri[i]][0]], sites[triangles.vertices[pts_tri[i]][1]], sites[triangles.vertices[pts_tri[i]][2]])\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[1 0 3]\n",
        " [2 4 0]\n",
        " [2 1 0]\n",
        " [2 4 3]\n",
        " [2 1 3]]\n",
        "[[ 0.98656691  0.33351093]\n",
        " [ 0.44928433  0.558523  ]\n",
        " [ 0.54835299  0.11849635]\n",
        " [ 0.20050494  0.71994957]\n",
        " [ 0.53645758  0.06919823]]\n",
        "(array([ 0.5,  0.5]), 1, array([2, 4, 0], dtype=int32), array([ 0.54835299,  0.11849635]), array([ 0.53645758,  0.06919823]), array([ 0.98656691,  0.33351093]))\n"
       ]
      }
     ],
     "prompt_number": 166
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Interpolate the proper value from the triangle vertices.\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 171
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}